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1.1  Calculation  of  Lunar  Phase  and  Illumination 


Initiator:  I.M.  Hussey 

Project  No.:  0001  Problem  No.:  4510 

Using  spherical  trigonometry,  the  elongation  angle  and  the  phase  of 
the  moon  can  be  calculated  from  the  given  values  of  the  declination  angles 
and  hour  angles  of  the  sun  and  the  moon. 


ZHs0S  = Solar  Declination  (b) 

4Hm0M  = Lunar  Declination  (c) 

✓ H OH  = Difference  of  Hour  Angles  of  the  Sun  and  the  Moon  (A) 
s m 

^•SOM  = Elongation  (e) 

In  the  spherical  triangle  SNM  with  the  centre  of  the  sphere  at  0, 
the  elongation  angle  is  given  by 

cos  e = sin  b sin  c + cos  b cos  c cos  A 
The  phase  of  the  moon  is  related  to  the  elongation  angle  by 

phase  = h (1-cos  e) 

For  a given  phase  of  the  moon,  or,  equivalently,  a given  elongation 

angle,  the  illumination  of  the  moon  can  be  calculated  by  means  of  a 

curve-fitting  method  using  a set  of  known  elongation-illumination 

data  |x.jj  as  published  in  "Explanatory  Supplement  to  the  Ephemeris"!^ 

The  method  used  for  curve-fitting  is  based  on  a cubic  equation  fit 

algorithm  which  provides  the  coefficients  A , B , C , D for  each 

range  (x.  , x.  ..)  of  elongation.  The  illumination  F(x)  is  then 
i i * 3 2 

given  by  F(x)  = A^x  + x + x + D^  where  x^<x  <x^  . These 
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calculated  coefficients  for  the  cubic  equation  fit  are  given  in  the 
following  table: 


Range  of  Elongation 

A 

B 

C 

D 

0 

to 

20 

-.9748  xlO3 

.9674xl02 

0.0 

0.0 

* 

20 

to 

40 

-.3425  xlO 

.1739x10 

.2810x10 

-.  2238x10" 1 

40 

to 

60 

-.1340 

-.3453 

.2992x10 

.1 674x1 0_1 

1 

60 

to 

80 

.3842  xlO-1 

-.1907 

.3015x10 

-. 6900x10" 1 

80 

to 

100 

-.1321  xlO-1 

.1652 

.2225x10 

.4851 

♦ 

100 

to 

120 

-.1227  xlO"2 

. 1926x10" 1 

.2761x10 

-. 3800x10" 1 

120 

to 

140 

.1193  xlO-1 

-.3982 

.7101x10 

-. 1487xl02 

140 

to 

160 

.2652  xlO"3 

.2003X10"1 

.2344x10 

.1717x10 

160 

to 

180 

-.7256  xlO"2 

.5625 

- . 1058xl02 

. 1036xl03 

Existing  solar-lunar  ephemeris  programs  (Project  0001;  Problems  1131, 
1461)  were  adapted,  and  output  topocentric  viewing  angles  as  well  as  the 
phase  and  illumination  at  specified  time  intervals. 


1.2  STELAZ  and  Polaris 
Initiator:  I.M.  Hussey 

Project  No.:  0001  Problem  No.:  4517 

Given  the  mean  place  of  a star  at  any  epoch,  e.g.  1975.0,  January  0 .978 
in  Right  Ascension  and  Declination  («nt0,^’0);  find  the  Apparent  Right 
Ascension  and  Declination  for  any  specified  date  and  time  by  correcting 
for  aberration  (primarily  stellar  annual),  precession,  and  nutation. 

For  any  geographic  location  (geodetic  latitude,  longitude,  and  altitude) 
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find  the  stellar  azimuth  and  zenith  angle  corrected  for  refraction. 

Obtain  culmination  and  elongation  angles  and  times  for  a selected  star. 

1.2.1  Functional  Description 

Following  the  "Explanatory  Supplement  to  the  Ephemeris"  P.  150^,  the 
above  three  corrections  are  applied  in  the  order  specified.  Annual 
parallax,  proper  motion,  and  orbital  motion  (of  star)  are  neglected. 

From  apparent  right  ascension  of  station  at  U.T.  hours  and  the  geodetic 
latitude  obtain  the  azimuth  and  unrefracted  elevation.  Correct  for 
refraction,  assuming  station  is  near  sea  level.  Upper  and  lower  culmin- 

0 O 

ation  angles  and  times  occur  for  azimuth  0 and  180  , but  not  necessarily 
respectively.  Linear  interpolation  is  used.  W and  E elongations  are 
determined  from  a parabolic  fit  to  3 points  near  the  extrema. 

NOTE:  The  spherical  geometry  calculations  show 
negligible  loss  in  accuracy  even  for  nearly  cir- 
cumpolar stars.  If  declination  >89.8  0 a warning 
message  is  printed. 

Figure  1 shows  a sample  computer  print-out  for  observation  of  four  stars 
at  4 hour  intervals.  Daily  culmination  and  elongation  angles  and  times 
have  been  requested  for  Polaris  and  are  also  printed.  Output  is  suppressed 
if  the  star  is  unobservable. 

The  program  has  been  verified  using  epochs  and  mean  places  of  stars  from 
different  annual  editions  of  the  American  Ephemeris  and  Nautical  Almanac, 
and  against  an  earlier  program  which  required  transcription  of  the  Pole 
Star  Table  from  these  Almanacs.  Results  are  reproducible  to  a few 
tenths  of  a minute  of  arc. 

1.2.2  Components 

Reference  Time  and  Coordinates: 

Dflh  Julian  day  at  0h  UT  = 2442731.5  (example) 

Tu  = (Doh-  2415020.  )/36525. 

TTu  = Tu  - .01  x AINT  (Tu  x 100.) 

u 

Mean  R.A.  of  Greenwich  at  0 UT 
<*G  = 6.645065  + 2400.  x TTu  + .05126  x Tu  (hrs.) 

= 24.065710  x (UT)/24 
DJ  = DQh  + UT/24. 
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Figure  1.  STELAZ  Output  with  Culmination  and  Elongation  Angles  for  Polaris 


u«poch 

-r 


= DQh  + UT/24.  - 2415020.0  Jan  Q5 

= d /36525. 

= Julian  day  number  for  1975.0  0 .978  (say) 
= (DJ  - D^poch)  / 365.2500 
= T - T/2. 


The  fundamental  argumentsA,  Vand  fcare  available  with  sufficient  accuracy 
from  Reference  3,  using  the  reference  times  determined  above. 

Aberration  Correction 

At  any  time,  find  Ethe  obliquity  of  the  ecliptic  and  v the  longitude 
of  the  sun  form  the  mean  equinox  of  date.  Constant  of  aberration 
K = 20.47".  Then,  jfrom  "Spherical  Astronomy"  by  W.M.  Smart,  p.  I84J 

»»_ 

C = -K  cos  £ cos  V 
D = -K  sin  y* 

cot=  1/15  coso<sec<f  ; d^  = 1/15  sin«*  sec^ 
c^=  tantcos/-  sinocsin^  ; d^  = cos«csin/ 
and  apparent  right  ascension  and  declination  is  obtained  from  true 
right  ascension  and  declination  by: 

= oc+  C c^t  D d^ 


- <f  + C Cj  + D d^- 


Annual  Precession  Correction 


m = 3s  .07234  + 0s  .00186  T 

n = 2o". 0468  - 0"  .0085  T 
= Is. 33646  - 0s  .00057  T 


where  T is  in  tropical  centuries  from  1900.0. 
Reduction  from  mean  to  apparent  places: 


= T(m  + n sin«ctan^) 

= T(n  cos ac) 

evaluate  at  mean  T over  the  period  T years  plus  frac. 


Nutation  Correction 


In  Longitude 


In  Obliquity 


<4^=  -17  .2327  * sin  (A) 


9 .21 


cos  (A) 
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where 


-A-  = Longitude  of  moon's  ascending  node 

= .71995354167  - .000147094228332  dr(l-. 2941759  x 10‘15dJ 

v L 

revolutions. 


where  d * Julian  day  number  - 2415020.0.. 
c 11900  .Sill 

£ o = (cost  + sintsin  ot  tan  ) Alf/-  cosoitan/At 

A& n s sintcosot  Av|>  + sine*  At 

1.2.3  Procedure  _ 

W (if) 

Apparent  (corrected)  R^A.  and  Declination: 


~1  a *.+  D d^ 

^ * So  + C^.  + D d^- 

«<p  =ot^  +T(m  + n sinotj  tan/ 

4p  = ^+  T (n  cosotj) 

o*  =o<p  + (cost+  sintsino<p  tan/p)^ij; 

- cosocp  tan/p^g 

<r-<fP  + sin£,  coscepA*j>+  since  p At 

Apparent  R.A.  of  station  at  UT(hrs.) 

©C  * 0.9174*  A +A*C+*r  - GL0NW  * (-^  ) 

S SeC  Q G (West)  dft9* 


Unrefracted  Azimuth  and  Zenith  Angle: 

Observer  station  orientation  on  the  celestial  sphere:  o (^(geodetic  latitude). 


Find  Zenith  Angle  from 
cos  b + sin  a sin  b cos  C 
Then,  sin  A = sin  C . 

& cos  A = cos  a - cos  b cos  c 

sin  b sin  c 

cos/ cos  * sin  C,  sin<f-  sln^cos  ZA) 
( sin/sin^  + cos/cos^cos  C) 
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1.3  Program  ECLPLOT 


Initiator:  J.  Aarons 

Project  No.:  4643  Problem  No.:  4844 

Program  ECLPLOT  was  designed  to  depict  the  local  eclipse  circum- 
stances for  any  solar  eclipse  and  any  observational  station,  at 
selected  time  intervals. 

1.3.1  Functional  description 

Program  ECLPLOT  uses  punched  output  of  program  ECLPT  which  includes 
exposure  fraction  EXP  for  every  time  step  as  well  as  semi -diameter 
of  the  sun  SDS,  semi-diameter  of  the  moon  SDM,  angular  distance 
ANGDIS  and  position  angle  POSANG  of  the  moon  relative  to  the  sun,  and 
position  of  the  sun's  axis  PSAXS.  (Project  0001,  Problem  3501) 

Based  on  this  information  program  ECLPLOT  produces  plots  which  depict 
the  progression  of  the  eclipse  and  also  indicates  north  and  sun-axis 
directions.  The  sun-moon  center  line  is  always  horizontal  to  allow 
a large  uniform  plot  scale  of  1 inch  for  5 secs  of  arc. 

1.3.2  Procedure 

Figure  2 shows  a typical  eclipse  configuration  and  identifies  the  variables 
described  above. 


Figure  2.  Solar  Eclipse  Viewing  Geometry 


Punched  card  output  and  format  is  as  follows: 


IYR,  IMO,  IDA,  IHR,  IMN,  EXP,  SDS,  SDM,  ANGDIS,  POSANG,  PSAXS 
Format  (5(12, IX) , 6F10.7) 

Exposure  EXP  is  a fraction  (1-obscuration),  Others,  SDS  to  PSAXS,  are 
in  radians.  Figure  3 shows  a typical  plotted  output. 


1.4  Program  SOLPLT 

Initiat-"-:  G.  Best 

Project  No.:  7635  Problem  No.:  4824 

In  planning  rocket  launches  or  other  experiments  depending  on  sunrise, 
sunset  or  twilight,  it  is  useful  to  have  plots  on  which  these  times 
are  given  for  a for  a specified  location.  These  plots 

can  be  made  even  more  informative  by  also  including  times  of  moonrise 
and  moonset. 

The  Solar  and  Lunar  Ephemeris  programs  previously  developed  and 
modified  under  project  0001,  Problem  Nos.  1131  and  1461  were  adapted 
for  this  problem.  Certain  refinements  also  make  this  SOLPLT  program 
more  up-to-date  than  the  original  programs. 

1.4.1  Functional  Description 
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5'  of  arc  * 1 inch 


EXP=0.5607  UT=  2010 

HAMILTON 

Figure  3.  ECLPLOT  Output  for  Solar  Eclipse  of  July  10,  1972 


Lj-  M 


Program  SOLPLT  determines  times  of  specific  elevation  and/or  depression 
angles  of  sun  and/or  moon  for  any  location.  A tabulation  and  a 
file  of  daily  values  is  created,  and  a plot  of  graphs  to  depict 
these  conditions  is  produced. 

(4) 

Program  SOLPLT  uses  program  SOLLUN  as  a principal  subroutine 
which  generates,  as  a function  of  the  universal  time,  the  hour  angles, 
azimuth,  declination  and  elevation  of  the  sun  and  the  moon  for  any 
station  location.  Program  SOLPLT  uses  the  above  subroutine  to  det- 
ermine local  standard  times  for  specified  elevation  and/or  depression 
angles  for  the  sun  and  the  moon.  Some  refinements  to  the  original 
programs  were  incorporated.  Double  precision  is  not  required  on  the 
CDC-6600  and  has  been  eliminated.  The  correction  for  ET-UT  is 
calculated  by  default  as  (KYR-1929)  secs.  All  parts  of  the  separate 
solar  and  lunar  programs  have  been  recombined  into  the  original  order 
(parts  A-M).  Refraction  correction  and  librations  of  the  moon  (part  M) 
are  not  used.  As  input,  program  SOLPLT  requires  only  the  coordinates 
of  the  station,  the  desired  values  of  the  elevation  and/or  depression 
angles  and  the  year  of  observations.  In  order  to  obtain  local  standard 
time  values,  the  applicable  standard  time  longitude  must  also  be  input. 
Note  that  for  angles  equal  to  zero  the  times  of  sun/mpon  rise  and  set 
will  be  produced. 

Subroutine  LORAN  uses  values  obtained  by  SOLPLT  to  produce  plots  of 
elevatio"  angle  times  which  are  represented  by  solid  lines  and  depres- 
sion angle  times  which  are  represented  by  broken  lines.  The  plot 
covers  a full  year.  Figure  4 shows  typical  plots  where  the  solar 
elevation  angle  contours  are  20°  apart. 

1.4.2  Mathematical  ur  Procedures 

To  calculate  the  time  for  the  specific  elevation  angle  the  following 
procedure  was  developed.  For  every  day  of  the  specified  year  the 
standard  time  monotonically  increases  starting  from  zero.  Using 
spherical  trigonometry  formulae  one  calculates  the  minimum  earth 
rotation  and  corresponding  increment  in  standard  time  to  step  from 
one  specified  angle  to  the  next.  Standard  time  should  increase 
monotonically  but  allowance  is  made  for  small  negative  increments 
to  avoid  possible  skipping  of  angles. 
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Figure  4.  Solar  Elevation  and  Moon  Rise-Set  Plots  for  1976  at  AFGL  and  20°  further  North 


For  elevation  angles  of  the  moon  the  program  uses  the  same  procedure 
as  for  the  sun,  with  subroutine  SOLLUN  providing  elevation  and  hour 
angles  of  the  moon. 

Whereas  rising  solar  elevations  always  occur  between  0 and  12  hours 
local  time,  and  dropping  elevations  after  12  hours,  this  property 
cannot  be  used  for  the  moon.  Instead,  a test  for  increasing  or 
decreasing  moon  elevation  is  made.  If  the  first  time  corresponds 
to  a decreasing  elevation,  the  moon  set  and  rise  times  for  the  day 
are  swapped. 
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Section  2.  Arctic  Ionospheric  Research 


2.1  Program  ISOZEN 


i 

Initiator:  J.  Buchau 


T 


Project  No.:  7663  Problem  No.:  4714 

ISOZEN  plots  contours  of  constant  solar  zenith  angle  for  any  90°  polar 
segment  of  the  Northern  Hemisphere.  The  contours  are  mapped  in  either 
geodetic  or  geomagnetic  coordinates,  and  are  designed  to  be  used  with 
an  AF6L  world-map  overlay.  Figure  1 and  2 are  examples  of  such  plots 
for  different  quadrants. 


2.1.1  Functional  description 


For  a given  date,  time,  quadrant  and  coordinate  system,  (geodetic  or 
geomagnetic),  ISOZEN  produces  a 10"  one-quadrant  plot  of  constant 
zenith  angle,  each  contour  being  separated  by  5°,  and  drawn  as  a 
solid  or  dashed  line,  solid  for  less  than  90°,  dashed  otherwise.  The 
present  version  allows  one  to  step  the  time  by  a given  increment, 
thereby  generating  several  plots  in  one  run. 

2.1.2  Mathematical  & Logical  Procedures 

The  solar  zenith  angle,  Z,  (the  angle  between  the  local  vertical  and 
the  position  of  the  sun),  is  a function  of  latitude,  longitude,  and  time. 

Z = Z (*,  $,  T) 

A contour  of  constant  zenith  angle  for  a particular  time  is  given  by  the 
equation 

Zo  = Z (a,  $,  To) 

where  Zo  and  To  are  given  constants. 

The  following  procedure  was  used  to  generate  such  contours » making  use  of 
some  subroutines  develoiW  previously'^ : 

A search  is  initiated  along  one  of  the  quandrant  boundaries  for  desired 
contour  values.  When  one  is  found,  either  6-  or  0 is  stepped  away  from  the 
boundary,  and  a Mueller’s  iteration  technique,  using  the  subroutine  SOtUN 
to  evaluate  Z,  is  used  to  find  the  other  coordinate  (0  or  0). 

This  procedure,  analogous  to  a blind  man’s  use  of  his  cane  to  navigate 
crosstown,  is  repeated  until  another  boundary  is  encountered,  thereby 
completeing  this  particular  contour.  The  intersection  is  catalogued  to 
avoid  redundant  plotting,  and  entire  sequence  repeated  until  all  contours 
have  been  mapped.  The  transformation  from  geodetic  to  geomagnetic  coord- 
inates, if  required,  is  effected  by  the  subroutine  CGINV. 
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Figure  2. 


2.2  Photometric  Data  Presentation 


Initiator:  J.  Buchau 

Project  No.:  5631  Problem  No.:  4605 

Program  PDTRK  is  essentially  a plotting  program  to  display  voltage  readings  from 
an  aircraft  mounted  360  degrees  scanning  photometer  along  the  projection  on  the 
spheroid  of  the  aircraft  flight  path,  all  in  a corrected  geomagnetic  polar  coordin 
ate  system. 

The  voltage  readings  are  plotted  around  an  approximate  circle  with  the  aircraft 
position  at  center,  the  circle  being  the  projection  on  the  spheroid  of  the 
intersection  of  the  scanning  ray  and  a shell  at  a given  altitude,  representing 
the  height  of  the  aurora. 

Early  versions  of  PDTRK  plotted  the  voltages  along  either  side  of  the  velocity 
vector  in  non  overlapping  semi  circles,  and  produced  similar  plots  looking  in 
the  opposite  direction.  These  were  found  insufficient,  however,  and  in  the 
final  version,  the  entire  circle  of  voltages  was  displayed  at  five  minute 
intervals,  along  with  a velocity  vector  from  the  aircraft  position  and  time/ 
date/ frequency  channel  information.  In  this  final  version,  the  background 
coordinate  system  was  restricted  to  the  portion  immediately  surrounding  the 
circle  of  values,  the  plots  were  centered  on  the  aircraft  position,  and  they 
were  all  put  on  35mm  film  strips  so  that  examining  the  plots  sequentially  could 
be  done  quickly.  Figures  3 and  4 show  two  successive  plots. 

2.2.1  Functional  Description 

Input  into  PDTRK  consisted  of  time  and  position  information  taken  from  the 
aircraft  log  and  entered  punched  card  data,  along  with  a magnetic  tape  containing 
the  photometer  voltages  for  the  flight,  which  were  read  at  approximately  every 
thirty-four  seconds.  The  aircraft  log  position  data  was  taken  at  variable  inter- 
vals which  ranged  from  ten  to  twenty  minutes  in  length. 

In  order  to  produce  the  plots  at  five  minute  intervals,  it  was  necessary  to 
interpolate  from  the  log  data  for  the  aircraft  position,  and  to  position  the  data 
tape  to  obtain  the  voltages  for  the  required  times.  The  data  from  the  closest 
time  was  used,  since  the  error  was  no  more  than  about  fifteen  seconds. 

Positioning  of  the  original  data  tape  proved  to  be  cumbersome,  so  a preprocessor 
program  PHOTO  was  written  to  produce  records  that  were  one  hundred  characters 
plus  8 character  fill  in  length,  ten  per  photometer  sampling. 

To  interpolate  the  aircraft  position  between  successive  log  points,  a subroutine 
POS  was  written.  Given  two  sets  of  time-positive  coordinates  and  an  intermediate 
time,  POS  will  return  the  position  of  the  corresponding  intermediate  point  on 
the  great  circle  between  the  given  points,  assuming  a constant  velocity  over  the 
entire  interval. 

To  find  the  projection  of  the  photometer  ray  and  the  aurora,  program  SAP  was  em- 
ployed. C0RRGM2  converted  geographic  coordinates  to  corrected  geomagnetic 
coordinates.  SAP  and  C0RRGM2  are  described  elsewhere'1' 

2.2.2  Mathematical  or  Logical  Procedures 


Subroutine  POS  is  essentially  the  solution  of  two  spherical  triangles. 


In  the  diagram  the  first  application  is  to  the  triangle  PjOPg.  where 
^ are  the  9^ven  points,  0 is  the  pole  and  OPj,  0P2  are  co-latitudes, 

to  obtain  P^  and  the  spherical  angle  OP^.  PjP  is  obtained  by  a 
simple  ratio  since  velocity  Is  constant,  and  then  the  spherical  triangle 
OPjP  Is  solved  to  obtain  OP,  the  new  co-latitude  and  PjOP,  the  longitude 
Increment. 

P 


To  produce  the  restricted  background  plot,  we  obtained  the  limits  of 
the  polar  coordinates  of  the  projected  number  field,  and  chose  conven- 
ient latitudes  and  longitudes  to  enclose  the  figure.  The  direction 
| j arrow  was  added,  and  the  rectangualr  coordinates  for  plotting  were 

corrected  by  the  coordinates  of  the  aircraft  position  to  center  the 
plot  on  the  aircraft. 


Section  2.3  DAASM  Flight  Program 
Initiatpr:  J.  Videberg 

Project  No.:  4141  Problem  No.:  4693,4748 

Introduction 


In  the  Blue  Deck  Polar  Model  version  of  the  ITS  78  Ionospheric  Predictions 
Program,  the  mode  7 output  tabulates  ionospheric  perameters  for  MUF  and  selected 
frequencies  for  a given  receiver  and  transmitter  configuration  repeating  at 
designated  steps  over  some  time  interval.  The  program  was  modified  so  that  in 
this  output  mode  aircraft  flight  log  information  could  be  entered,  effectively 
causing  transmitter  coordinates  to  vary  as  well  as  the  time  coordinate.  Then 
using  interpolation  methods  the  ionospheric  perameters  are  output  for  selected 
equal  time  intervals  over  the  aircraft  flight. 

In  addition  a table  is  presented  between  each  two  time  points  which  summarizes 
the  Doppler  frequency  shift  for  each  mode  and  each  frequency. 

Functional  Description 

The  calls  to  subroutine  CGMTIME  were  all  set  up  for  times  which  were  integral 
hours.  Since  we  wished  to  handle  random  times,  these  calls  had  to  be  reworked 
to  permit  hours,  minutes  and  seconds.  The  change  to  the  program  essentially 
involved  a)  eliminating  the  automatic  hourly  integral  time  increments  b)  per- 
mitting variable  transmitter  coordinates  c)  interpolating  d)  summarizing  Doppler. 

The  reorganization  of  the  program  was  accomplished  by  writing  a new  main  program, 
called  NUMAN  which  reads  input  information  and  reorganizes  it  into  files  called 
TAPE2  and  TAPE3.  With  the  original  time  steps  inhibited,  it  then  calls  the 
original  BDPM  main  program,  as  a subroutine. 

NUMAN  also  does  the  interpolation  along  the  flight  path,  utilizing  subroutine  POS, 
which  is  discussed  in  sec.  2.2  of  this  report. 

Upon  completion  of  outputing  the  ionosphere  parameters,  BDPM  calls  subroutine 
OUTS,  which  outputs  the  Doppler  summary  information.  The  printing  is  done  by 
subroutine  SUBPR,  which  is  called  by  OUTS. 

INPUT 

card  input 


cd 

1 

Time  increment 

Doppler  Summary 

cc 

Format 

Contents 

1-5 

15 

Time 

6-10 

15 

Doppler  summary  selector 

0- -no  summary 

1- -summary 

cd 

2,  3 

Transmitter  and  receiver  antenna  cards,  as  described  in  ESSA  TECHNICAL 
REPORT  ERL  1 10- ITS  78  (denoted  henceforth  as  ETR) 

cd 

4 

99  in  cc  7-8. 

cd 

5 

Month  and  Sunspot  card  (see  ETR) 

32 


r 


1 


cd  6 -1  In  cc  4-5 

cd  7 1 Program  control  cards  (see  ETR) 

cd  7+k  f 

cd  8+k  99999  in  cc  1-5 


cd  9+k  Frequency  complement  card  (see  ETR) 
cd  10+k  Circuit  card  (see  ETR) 


cd  11+k 

iFlight  time-position  cards 

cc 

format 

contents 

1-5 

15 

Universal  time-hrs( integer) 

5-10 

F5.0 

Universal  time-min(decimal ) 

cd  11+k+l 

11-15 

15 

Latitude-deg. 

(integer) 

16-20 

15 

Latitude-mi n 

(integer) 

21-25 

15 

Longitude-deg. 

(integer) 

26-30 

15 

Longitude-mi  rj 

(integer) 

cd  12+k+l 

1000  in 

cc  2-5. 

Tape  input 

1) .  Tables  for  geographic-geomagnetic  coordinate  conversion. 

2) .  Numerical  coefficients  for  various  ionospheric  layers  as  described  in  ETR. 
OUTPUT 

1).  Mode  7 output  as  described  in  ETR  with  varying  transmitter  coordinates. 


2).  Doppler  rrmmary  table  (see  Fig.  5).  For  each  mode-frequency  box,  in 

printed  beginning  time  delay,  ending  time  delay,  difference,  and  Doppler 
frequency  shift. 

JOB  CONTROL 

Organizing  the  modified  program  so  that  the  imput  could  be  read  as  called  for  by 
the  subroutine  was  accomplished  by  job  control  cards. 

The  coordinate  conversion  tables  were  on  a mass  storage  permanent  file,  as  were 
the  coefficients.  With  the  former  designated  BDAT,  and  the  latter  BDPDAf  and 
the  binary  program  file  as  BPR,  the  following  job  control  was  used. 

ATTACH,  BB,  BPR,  ID=  — . 

ATTACH,  TAPE  9,  BDPDA,  ID=  — . 

ATTACH,  CC,  BDAT,  ID=  — . 

COPYCR,  CC,  EE. 

COPY  CR,  INPUT,  EE. 

REWIND,  EE. 

COMBINE,  EE,  FF,  2. 

BB,  FF,  PL =77777. 
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Figure  5.  Doppler  Summary  Table  for  5 Minutes  of  Flight  Time 


Initiator:  C.  Pike 
Project  No.:  7663 


2.4  ISIS-2  Data  Base 


Problem  No.:  4604 


Three  magnetic  tapes  containing  ISIS-2  satellite  soundings  of  electron  density 
profile  versus  height  were  received  from  the  Canadian  Communications  Research 
Center  and  were  reduced,  sorted  and  merged  into  one  data  base  tape  for  sub- 
sequent research  at  AFGL.  The  tapes  were  converted  from  9 track,  800  bpi 
EBCDIC  format  to  SCOPE  7 track  800  bpi,  BCD  format.  The  source  tapes  were 
in  the  "loge  (density)"  format  that  was  also  used  for  the  later  Alouette  I N( h 1 
data.  Program  FIRSTFL  had  previously  been  developed  for  processing  that  data(D 
and  was  modified  to  generate  the  ISIS-2  data  base,  including  geographic,  geomag- 
netic, and  solar  environmental  information.  Over  25,000  samples  were  obtained 
after  rejection  of  between  1 and  2 percent  because  of  various  unacceptable  data 
conditions.  The  electron  density  is  provided  at  altitudes  from  950  km  to  350  km 
in  steps  of  100  km.  For  altitudes  below  Hmax  the  electron  density  is  entered  as 
-9.99.  Maximum  electron  density  EDmax  is  replaced  by  fgF2  according  to  EDmax= 
1.24*104(foF2)2.  The  data  were  time  sorted  except  where  the  F5.2  format  for  UT 
did  not  provide  enough  resolution.  The  data  covers  1970  thru  1973.  Table  1 gives 
the  format  for  the  tape  that  was  created. 
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Table  1. 


ISIS-2  Data 


i 

i 

| 


Variable  Description 
ID  (data  source=ll) 

Year 

Day  Number  (1-366) 

UT  (hours) 

Geog.  Latitude 
Local  Time  (hours) 

Mag.  Latitude 
Mag.  Time  (hours) 

Geog.  Longitude-  E positive 
Cosine  of  Zenith  Angle 
Kp' 

F10.7cm  Solar  Flux 

Quality  figure  (4- very  good,  9- very  poor) 
Altitude  of  Satellite  (km) 

Elec.  Dens,  at  Satell itexlO"^ 

TOTALN  x 10-11 

Elec.  Dens.  (950,  850.... 350)  x 10"^ 
f0F2 

^max^m) 


Format 

12 

12 

13 

F5.2 

F6.2 

F5.2 

F6.2 

F5.2 

F5.1 

F6.3 

F4.2 

13 
11 

14 

F6.2 

F8.3 

7F6.2 

F4.1 

F4.0 
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Initiator: 
Project  No 


2.5  Auroral  Morphology  Study  (Part  II) 


R.  Nadile 

. : 7670  Problem  No.:  4689 


Program  AIMS,  which  provides  a CRT  plot  satellite  view  of  the  aurora  and  also 
depicts  boundaries,  cross-sections,  and  solar  aspect,  (1)  was  modified  for 
extended  uses. 

1. )  A polar  projection  of  the  earth  and  auroral  oval  including  a map  of 
continental  outlines  was  programmed.  Intersection  of  aurora  and  satellite  field 
of  view  is  depicted,  with  contributions  to  a selected  tangent  height  also  shown. 

2. )  An  improved  sunline  routine  was  developed  to  draw  the  terminator  on 
these  plots. 

3. )  The  program  was  adjusted  to  also  allow  southern  hemisphere  auroral  morphology 
studies. 

2.5.1  Polar  Projection 

Figure  6 shows  a typical  CRT  plot  for  a geo-stationary  satellite  with  emission 
profile  betweeen  75  and  150  km  as  a function  of  altitude  h given  by: 

V = HN  x e ( 1-HNl* !6) 


where 


HN  = (h  -75  ) / 15 . 

The  auroral  oval  is  taken  to  extend  uniformily  betweeen  61.5°  and  73°  N;  geo- 
magnetic latitude. 

The  polar  piojection  is  designed  to  show  the  viewing  geometry,  subauroral  coordin- 
ates, and  intensity  of  contribution  to  the  visible  aurora  for  a specified  tangent 
height.  Continental  outlines  may  be  included.  Figure  7 shows  such  a CRT  plot 
for  a 30  km  tangent  height  section  for  the  example  of  Figure  6.  Since  30  km 
above  the  earth  limb  is  the  lowest  tangent  height  of  interest,  the  shading  at 
higher  elevations  would  progressively  cover  the  two  small  regions  bounded  by  the 
aurora.  The  rest  of  the  aurora  remains  invisible,  either  in  front  of  or  behind 
the  earth  limb. 

The  limits  of  azimuth  AZ  and  elevation  EL  are  converted  to  a dashed  viewing  line 
on  the  paper  from  the  satellite  located  on  the  left.  Assuming  the  satellite  is 
at  latitude  0,  the  slope  of  the  corresponding  dashed  line  is  given  by 

ci  nnc  tan  EL  sin  AZ 

iLUKt  ' ‘ COS  0 + tan  EL  COS  AZ  sin  0 


2.5.2  Solar  Terminator  at  Shell  Layer  Height 

The  sunline  generating  routine  developed  previously^  moved  the  terminator  on 
the  earth's  surface  back  by  a distance  SDL  from  the  direction  of  the  solar 
radius  vector  § j,  5,.  In  this  satellite-oriented  geocentric  coordinate 
system,  the  x-axis  is  to  the  north,  y to  the  east,  and  z passes  through  the  satel- 
lite. Equal  increments  were  taken  in  the  x-direction  and  the  terminator  was 
solved  for.  Difficulties  were  occasionally  encountered  when  the  sun  vector  was 
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close  to  the  x-z  plane.  An  improved  routine  has  been  implemented 
in  which  the  farthest  point  of  the  terminator  on  and  behind  the  earth 
is  first  determined,  and  then  the  vector  to  this  point  is  rotated 
about  the  earth  sun  axis  in  * D ' angle  steps.  The  sign  of  the  rotation 
depends  upon  the  sign  of  f,  . 

Thus,  starting  with 

5 +‘ 

solve  for  5, as  before. 

Then,  for  each  step 


(j.i  -i.i ) 

(■U  • 5,5 ) 

-yi) 


0 

cos  D 
-sin  D 


This  terminator  is  moved  back  on  to  the  shell  by  the  shadow  distance  SDL. 

S5  = K-  SDL  Y)s=rl-  SDL  2.,  §s  = SDL  . 

Finally,  if  ?.s<5jthe  terminator  on  the  shell  is  obscured  behind  the 
earth  limb. 

is  given  by  ^ - SDl/P2-1 


where  $ - earth-satellite  distance  in  earth  radii. 


§ 

f 


2.6  Ionospheric  Research  Program  Library  UPDATE  Tape 


A group  of  programs  that  were  developed  over  a number  of  years  in 
support  of  ionospheric  research,  and  that  have  some  commonality  either 
as  regards  use  or  subroutines  were  consolidated  into  an  UPDATE  program 
library.  Operation  consists  of  having  a binary  subroutine  library 
available,  created  under  EDITLIB,  and  compiling  and  executing  only 
the  requisite  main  program.  The  following  programs  comprise  the  present 
library: 


PROGRAM 


1)  TRK2 

}2)  ALTRAK2 
3)  STATS 
4)  CON 
5)  CGARC1 
6)  AIMS 
7)  A I MSA 
8)  CONFMAP 
9)  WRLDMAP 
10)  MJRECT 
11)  ABCD 
12)  ISOZEN 


Flight  Tracks 

Satellite  Tracks 

Statistics 

Contour  plots 

Auroral  are  plots 

Auroral  Morphology  (Earth  Limb) 

Auroral  Morphology  (Global) 

DAPP  - CGM  grid 
Continental  Outlines 
Rectangular  Coord.  CGM  plot 
Create  Binary  CGM  file 
Iso-zenith  angle  contours 


These  progr-ms  are  described  in  an  earlier  report'  'or  in  the  present  one. 


Reference 

1)  "Analysis  and  Programming  for  Scientific  Research",  Logicon,  Inc.,  Final 
Report,  AFCRL-TR-74-0480,  September  1974. 


3.0  Rapid  Orbit  Prediction  Program  - ROPP 


Initiator:  I.  M.  Hussey 

Project  No.:  0001  Problem  No.:  4506,  4801 

For  drag  perturbed  orbits  analytic  perturbation  theory  excluding  short  period 
oscillations  is  used  in  program  ROPP  with  large  step  sizes  to  provide  reasonably 
accurate  and  rapid  prediction  of  the  orbit  elements. ^ This  program  was 
modified  to  include  additional  features: 

1)  Estimation  of  the  ballistic  coefficient  CpA  relative  to  the  atmospheric  density 

U 

model  based  on  two  sets  of  elements. 

2)  Input  and  output  of  the  osculating  position  and  velocity  vector. 

3)  Altitude,  perigee-apogee,  eclipse  in-out  times,  revolution  number 
information,  and  apogee-perigee  plot  file. 

The  present  program  typically  outputs  parameters  1-4  times  per  second  on  the 
CDC6600.  For  example,  a low  altitude  satellite  run  with  print-out  every  12  hours 
executes  1-2  days'  of  orbit  per  second. 


3.1  CpA  Optimization 
im 

Program  ROPP  was  designed  to  provide  an  approximate  satellite  ephemeris  and  life- 
time, given  elements  and  drag  factor  at  any  epoch.  For  the  drag  model  the 
geomagnetic  anc  solar  flux  activity  are  required  as  well  as  the  ballistic  coefficient 
CqA/m  where  CQ  is  the  drag  coefficient  and  A/M  is  the  area  to  mass  ratio.  For 

satellites  at  significant-drag  altitudes  the  mean  motion  and  resulting  in-track 
error  is  extremely  sensitive  to  inaccuracies  in  the  atmospheric  density  model 
or,  what  is  essentially  the  --.e  thing,  the  ballistic  coefficent.  Given  two 
sets  of  elements  plus  revolution  number  reasonably  separated  in  time,  CDA/M 

may  be  adjusted  in  ROPP  to  cause  the  true  anomaly  to  correspond  closely  to  that 
of  the  second  element  set.  The  semi-major  axis  (A),  eccentricity  (E),  inclination 
(I),  argument  of  perigee  (W  or  0),  and  ascending  node  (N)  are  slowly  changing 
elements,  so  that  great  confidence  can  then  be  placed  in  the  trajectory  between 
the  element  times  and  for  a considerable  period  beyond. 

The  main  problem  in  implementing  this  optimization  procedure  is  in  computing  the 
cumulative  angular  rotation  from  the  first  element  set  to  the  second.  After  a 
few  weeks  the  average  mean  motion  is  unable  to  predict  the  right  revolution  number. 
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It  is  thus  necessary  at  each  incremental  output  time  to  correct  the  cumulative 
angular  rotation  obtained  from  the  mean  motion  and  rate  of  change  of  argument  of 
perigee,  so  as  to  match  the  true  anomaly  obtained  from  the  elements  at  that  time. 
This  integration  can  then  be  safely  continued  up  to  the  second  set  of  elements. 
The  procedure  is  described  below: 

Let  I be  any  intermediate  time,  and  X the  next  incremental  time.  Let 
1 and  2 denote  the  initial  and  final  element  sets  and  times. 


Notation: 


N - rev.  number 


v - true  anomaly 
w - argument  of  perigee 

4 

w - rate  of  change  of  w 
MM  - mean  motion 

AT  - cumulative  rotation  from  equator  crossing  preceding  initial  elements 
by  (vl  + wl). 

a)  If  ATX  is  the  cumulative  rotation  at  time  X the  corresponding  revolution 

number  is  given  by  NX  = NA  + I NT  [(ATX/ 360  f]  (1) 

Where  NA  = N1  - INT  £(vl  + wl)/36oj  (2) 

NA  gives  the  base  revolution  number  to  be  added  to  any  subsequent  calculation 
to  give  the  correct  rev.  number. 

b)  Compute  rotation  increment  for  interval  ^t  from  I to  X 
Estimated  =^MMI+MMX^ ^t  + |wI+wX^  At 

Calculate  ATX  = ATI+AA 

Adjust  ATX  to  match  (vl+wl)  to  the  closest  360*. 

c)  This  ATX  becomes  ATI  for  the  next  step;  continue  for  each  increment  t 
until  time  2 where  we  get  AT2,  the  predicted  revolution  number 

NT  = NA  + INT  (AT2/360)  (3) 

and  the  predicted  total  rotation AAT  = AT2  - ATI.  (4) 

d)  Calculate  actual  rotation  DAT  from  time  1 to  time  2. 

DAT  = ( N2-N1 )x  360  + MOD  (v2+w2,360)-  MOD  (Vj+Wj.360)  (5) 

Finally,  error  between  observed  and  computed  rotation  is 

Error  = DAT-^AT  (6) 


jl 
! . 


If  the  error  Is  positive  C^A  must  be  increased.  Successive  iterative  experiments 

1 T 

with  CqA  adjusted  also  provide  the  rotational  change  sensitivity  to  these  ad- 

TT 

justments,  thus  allowing  a rapid  convergence  to  the  optimum  ballistic  coefficient. 
The  procedure  was  implemented  with  subroutine  CCDAM,  and  modification  of  existing 
routines  including  OUTSUB  and  READIN. 

Figure  1 shows  a typical  computer  output  where  CQA  was  optimized.  The  first 

IT 

element  set  is  entered  in  the  normal  way  and  is  now  shown.  The  NAMELIST  input  is 
(2) 

also  standard'  , except  that  the  optimization  procedure  is  initiated  by  setting 
both  CD  and  CDAM  to  zero.  The  second  element  set  is  shown  entered  as  a position- 
velocity  vector  after  the  $END  card.  Iteration  0 provides  a crude  first  estimate 
for  CDAM  based  on  a mean  motion  that  has  increased  from  time  1 to  time  2.  Three 
subsequent  trajectory  iterations  are  shown  before  a close  enough  match  is  obtained 
to  the  nodal  period.  The  final  ROPP  prediction  run  with  output  is  then  made  with 
CDAM*. 04247  ft.2/lb. 


3.2  Position  and  Velocity  I; 


The  basis  for  the  addition  of  osculating  position  and  velocity  input  - output 
was  the  gravitational  perturbution  model  that  exists  in  ROPP.  Subroutine  OSCSET 
generates  comprehensive  functions  of  eccentricity  and  inclination  plus  their 
derivatives.  These  are  used  by  subroutine  OSCSUB  to  compute  osculating  Kepler 
elements.  Corw^sion  between  these  elements  and  position-velocity  is  a standard 


procedure^  and  was  implemented  with  subroutines  ELMTPV  and  PVTELM^V 


Conversion 


from  osculating  to  mean  Kepler  elements  was  added  into  PVTELM  as  an  iterative 
procedure  where  the  mean  elements  are  corrected  by  the  difference  between  the 
desired  and  resulting  osculating  elements.  Four  iterations  including  OSCSET 
recalls  are  sufficient  foi  oarce.  The  mean  Kepler  elements  in  Figure  1 
were  derived  in  this  manner  from  SCF  position-velocity  elements. 

Figure  2 shows  the  normal  full  output,  in  this  case  twice  a day.  The  last  two 
lines  present  osculating  position  and  velocity  which  have  been  made  consistent 
with  the  original  input  due  to  the  procedure  described  above. 


3.3  Other  Modifications 

a)  Program  ROPP  was  converted  from  double  precision  (RUN)  to  single  precision 
(FTN)  with  negligible  loss  of  accuracy.  Core  memory  requirement  was  reduced 
by  30K  octal,  compilation  time  by  90  percent,  and  execution  time  by  about 
50  percent. 


Figure  1.  ROPP  CpA/M  Optimization  for  Elements  40  Days  Apart 
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b)  Element  input  may  be  ADC  2-card  or  5-card  sets  (see  Section  4.1)  or  SCF 
position-velocity  data.  These  capabilities  parallel  program  LOKANGL. 


c)  When  using  CQA  optimization,  if  the  desired  time  span  for  the  trajectory 

ir- 
is not  input,  the  program  operates  between  the  times  of  the  two  sets  of  elements. 


d)  Mean  satellite  altitude  is  included  with  latitude  and  longitude  output. 


If  R is  the  mean  radius  vector  and  0 the  latitude  of  the  satellite, 
the  altitude  of  the  satellite  H$  is  given  by 


a<y 


iLf^n  Sin  zt 

d-f)2 


where  afi  is  the  equatorial  radius  of  the  earth,  and  f is  the  earth 
flattening. 


e)  Local  times  at  perigee  and  apogee  are  included  in  the  output,  and  are 
readily  derived  from  existing  quantities. 


f)  Universal  times  of  entry  into  and  exit  from  solar  eclipse  are  included 
along  with  the  eclipse  duration,  for  a period  close  to  the  requested  output 
time. 

g)  Revolution  number  at  any  time  is  calculated  according  to  Equation  (1) 
of  Section  3.1,  and  is  output. 


h)  Optional  Apogee-Perigee  Plot  File 

If  a ' 1'  is  punched  in  column  71  of  the  first  card  of  the  first  element 
set,  apogee-perigee  information  is  written  on  TAPE7  in  punched  card  format  com- 
patible with  program  P LOUT  (see  Section  5.2).  This  output  is  described  in 
Table  1.  Figure  2 in  this  Section  is  part  of  the  run  for  satellite  8468  which 
subsequently  produced  the  apogee-perigee  plots  of  Figure  1 in  Section  5. 

In  another  application,  mean  elements  from  ROPP  were  punched  out  period- 
ically to  study  disparities  with  mean  elements  provided  by  ADC.  Results  of  this 
study  are  discussed  in  Section  4.3. 


3.4  Problems  with  ADC  and  SCF  Elements 

A number  of  problems  were  identified  during  the  use  of  the  orbit  prediction  pro- 
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Table  1.  ROPP  Apogee-Perigee  Output  Description 


VAR.# 

VAR. 

UNITS 

FORMAT 

CARD  COL. 

1 

Day  No. 

13 

1-3 

2 

HOUR 

12 

4-5 

3 

MIN 

12 

6-7 

4 

YR 

12 

8-9 

5 

LAT  (per.) 

Deg. 

F5.1 

10-14 

6 

WLONG(per. ) 

Deg. 

F5.1 

15-19 

7 

ALT  (per.) 

Km. 

F7.2 

20-26 

8 

LAT  (apog.) 

Deg. 

F9.1 

27-35 

9 

WLONG(apog-) 

Deg. 

F9.1 

36-44 

10 

ALT  (apog.) 

Km. 

F7.2 

45-51 

11 

REV.  No. 

19 

52-60 

12 

L.T.  (per) 

Hr. 

F6.2 

61-66 

13 

L.T.(apog.) 

Hr. 

F4.1 

67-70 

14 

MODEL 

"18" 

12 

71-72 

15 

14 

73-76 

16 

IDVEH 

14 

77-80 

! 

I 


grams.  Some  of  these  have  been  resolved,  and  others  that  would  require  further 
extensive  study  but  were  tolerable  are  described. 

Figure  3 shows  a printer-plot  of  the  differences  in  eccentricity  between  the  mean 

values  provided  by  ADC  (A),  and  the  mean  values  derived  by  conversion  of  the 

SCF  position-velocity  vector  (X).  The  element  sets  are  for  satellite  7499  and 

cover  a period  of  60  days.  The  argument  of  perigee  is  drawn  and  the  greatest 

o 

eA-ex  values  are  seen  to  occur  for  w=270  . A zero  drag  run  of  ROPP  showed  that 
this  program  inputs  and  outputs  mean  elements  that  include  the  long  period 
oscillation  due  to  J^.  The  larger  e^  value  (~*00075)  for  perigee  near  the  South 

pole  shows  that  the  ADC  elements  leave  out  the  long  period  perturbations  whereas 
ROPP  assumes  that  they  are  included.  The  ADC  elements  are  therefore  adjusted  for 
eccentricity  (see  Section  4.1,  Equation  1).  The  mean  semi-major  axis  is  obtained 
from  the  mean  motion  using  Equation  5 of  Section  13  of  Section  4.4.1. 

Figure  4 shows  the  same  mean  eccentricity  values  on  an  absolute  scale.  A ROPP 
pure  prediction  run  was  made  starting  with  a pair  of  early  elements  and  an  optimized 
C A 

D . The  mean  eccentricity  (R)  is  seen  to  drop  off  too  rapidly  after  30  days, 

M 

along  with  the  associated  perigee  altitude,  indicating  that  the  atmospheric  density 
model  in  ROPP  has  an  insufficient  scale  height,  (also  see  Figure  6,  Section  4). 

This  problem  was  not  tackled. 

Figure  5 represents  the  results  of  a search  for  the  cause  of  erratic  altitude 
predictions  in  ROPP  when  using  different  SCF  vector  sets.  The  reference  eccentricities 
(X)  are  the  ones  obtained  by  conversion  of  SCF  elements  in  LOKANGL  (Section  4.3). 

These  elements  were  also  converted  in  ROPP  and  the  differences  eR-ex  are  shown. 

The  ADC  mean  eccentricity  differences  eA~ex  are  also  shown.  The  abscissa  is  the 

argument  of  latitude  (w+M).  Values  for  the  different  element  sources  were  time 
interpolated  before  comparison.  The  ROPP  eccentricities  do  not  even  monotonically 
decrease.  The  conclusion  is  that  ROPP  conversion  between  mean  and  osculating 
elements  is  always  poor  when  these  elements  are  provided  near  the  poles.  Data  is 
absent  near  the  South  Pole,  but  a repeating  pattern  seems  indicated.  The  equally 
likely  large  positive  or  negative  differences--leading  to  the  erratic  altitudes 
is  unexplained.  Inspection  of  ADC  eccentricities  with  the  correction  also 
shows  substantial  disparity  with  the  SCF  derived  values;  however  the  fact  that 
ADC  elements  are  invariably  provided  near  an  upward  equator  crossing  limits  our 
evaluation  of  this  problem. 
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Figure  4.  ADC,  SCF  and  ROPP  (predicted)  Eccentricities 


Figure  5.  Comparison  of  ROPP  and  LOKANGL(X)  Derived  Eccentricities  Against  Argument  of  l.atitude 
ADC  eccentricities  are  also  compared. 


3.5  Parameter  - Lifetime  Study 

The  effect  of  varying  four  different  factors  independently  on  the  lifetime  prediction 
of  the  ROPP  program  was  studied.  The  same  position  velocity  vector  set  was  used 
in  every  case.  This  set  was  dated  25  days  before  the  observed  crash  date  for 
satellite  7499  (S3-1).  The  factors  and  their  effects  are  described  below.  Nominal 
values  for  these  parameters  were  CDAM=.055,  AP=15,  FL0BAR=75,  and  HMAX=.5. 

CDAM 

CDAM  is  the  least  known  of  the  four  parameters  studies.  Unfortunately,  it 
has  the  largest  effect  on  the  lifetime  prediction.  By  varying  CDAM  from  .02  to 
.08  one  finds  that  the  lifetime  prediction  varies  over  a range  of  45  days.  The 
ROPP  program  has  a sensitivity  of  about  7 days  for  every  .01  change  in  CDAM. 

AP 

The  AP  values  are  generally  well  known  and  a reasonable  average  value  over 
an  appropriate  time  span  is  used.  The  ROPP  program  is  sensitive  by  one  hour  for 
every  unit  change  in  AP.  This  effect  is  very  little. 

FlOBAR 

The  FlOBAR  values  are  also  well  known  and  are  used  just  as  the  AP  values. 

The  ROPP  program  lifetime  prediction  changes  by  approximately  one  hour  for  every 
unit  change  in  FLOBAR. 

HMAX 

HMAX  is  the  maximum  value  of  the  integration  step  size.  The  variation  of 
HMAX  from  .0625  to  1.0  produces  a change  in  the  lifetime  prediction  of  two  hours 
at  most.  The  only  important  consideration  in  choosing  a value  for  HMAX  is  that 
the  value  must  be  small  enough  to  prevent  large  errors  from  occurring  in  the 
calculations.  Generally  these  errors  show  up  as  a negative  perigee  altitude 
in  the  ROPP  output. 
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4.0  Satellite  Ephemeris  Generation  Program  - LOKANGL 
Initiator:  I.M.  Hussey 

Project  No.:  0001  Problem  No.:  4506,4801 


Program  LOKANGL  is  a specialized  orbit  determination  program  that  primarily 
serves  the  requirements  of  the  AFGL  research  community.  Radar  and  other  obser- 
vational data  may  be  dispensed  with  when  good  quality  orbital  elements  are 
provided  by  other  agencies  that  regularly  operate  full  scale  orbit  determin- 
ation programs.  The  program  analytically  generates  satellite  ephemerides  along 
with  related  observational  conditions  as  desired  by  various  researchers.  The 
program  has  been  substantially  augmented  in  response  to  growing  use.  LOKANGL 
was  modified  to  provide  reliable  ephemerides  over  long  spans  of  time  by  inter- 
polating between  sequences  of  elements.  Abnormal  discontinuities  due  to  improper 
elements  can  be  recognized  by  inspection.  The  program  has  also  been  extended 
to  accept  and  output  position-velocity  vectors  such  as  are  provided  by  SCF. 
Support  programs  for  plotting  key  orbital  characteristics  over  the  time  spans 
have  been  developed  and  are  regularly  used.  These  are  described  in  Section  5. 

4.1  Use  of  Element  - Sets 

Figures  1 and  2 describe  the  card  format  input  for  elements  that  are  accepted 
by  LOKANGL  and  R0PP.  Information  which  may  be  punched  but  is  ignored  by 
LOKANGL  is  shown  double-lined  out.  The  crossed-out  columns  are  blank  separators 
between  fields.  Figure  1 shows  the  original  ADC  format  which  was  received  by 
teletype.  The  sixth  line  was  not  provided.  The  one  digit  year  is  now  used 
in  LOKANGL  tc  indicate  a year  in  1970.  The  Modified  Julian  Days  on  card  2 is 
the  Julian  Date  plus  fraction  less  2400000.5  Thus  the  fraction  is  the  UT  in 
days.  The  integer  part  is  optional  and  is  shown  dashed.  In  early  1974  ADC 
switched  over  to  a 2-card  teletype  message  and  eliminated  the  semi-major  axis 
and  all  derivatives  ot  t.‘.v  ^'oLcn+s.  Only  the  first  and  second  derivatives 

• 9 0 

of  the  mean  motion  (N/2  and  N/6)  remain.  This  format  is  shown  as  Form  No.  2 
in  Figure  2.  A technical  document  ^ outlined. how  the  missing  quantities 
could  be  recovered.  This  brought  out  a problem  in  the  determination  of  the 
semi-major  axis  which  was  subsequently  resolved  (see  Section  4.4). 

Form  No.  3 in  Figure  2 is  adapted  for  entering  SCF  position-velocity  vector 
sets  that  are  also  available  by  teletype.  Period  Decay  in  seconds/rev.  provides 
the  drag  estimate.  Conversion  of  this  data  to  mean  Kepler  elements  for  use  in 
LOKANGL  is  described  in  Section  4.3. 
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Figure  1.  ADC  5-card  Element-set  Format 

Crossed  and  double-lined  entries  are  not  required 
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Figure  2.  ADC  2-card  and  SCF  Position-velocity  Element-set  Formats 


When  any  of  the  three  element  sets  are  used  in  ROPP,  derivative  information 
is  not  required  since  the  program  includes  a complete  force  model  and  analyti- 
cally generates  the  trajectory.  Program  LOKANGL  on  the  other  hand  operates 
primarily  by  polynomial  extrapolation  and  interpolation  of  mean  elements.  In 
the  case  of  SCF  data,  estimates  of  the  rates  of  change  of  the  elements  must 
be  made  initially  by  semi-analytic  means.  This  is  also  described  in  Section  4.3. 

When  successive  sets  of  elements  are  available,  rates  of  change  can  be  tested 
against  the  new  elements  and  adjustments  made.  This  is  discussed  in  Section  4.2. 


The  derivatives  provided  above  for  use  by  LOKANGL  cannot  include  the  long  period 

1 31 

perturbation  in  eccentricity  ' ' due  to  the  progression  of  the  argument  of 
perigee  over  the  northern  and  southern  hemispheres.  The  mean  mean  eccentricity 

is  therefore  converted  to  mean  eccentricity  by 

. _ - . .001066  . ...  ... 

e = e + sin  w sim  (1) 


before  computation  of  osculating  elements  in  LOKANGL.  In  the  case  of  SCF 
data  the  mean  eccentricity  is  initially  further  reduced  to  mean  mean  eccentri- 
city by  inverse  application  of  the  above  equation  6.  This  finally  makes  the  mean 
elements  from  the  three  sources  consistent  as  far  as  possible.  Some  results 
of  a study  that  was  made  to  resolve  these  discrepancies  are  discussed  in 
Section  4.4. 

4.2  Use  of  Successive  Sets  of  Elements 

Satellites  7499  (S3-1)  and  8468(S3-2)  have  been  of  particular  interest  because 
of  a number  of  AFGL  experiments  on  board.  ADC  and  SCF  element  sets  based 
on  their  full  scale  orbit  determination  programs  are  available  by  teletype  in 
a timely  fashion,  allowing  generation  of  ephemerides  using  a relatively  high 
density  of  elements.  In  fact,  since  the  telemetered  data  is  not  received  for 
some  time,  the  latest  elements  extend  well  beyond  the  available  data.  Program 
LOKANGL  was  modified  to  accept  time-successive  sets  of  elements  in  any  mix,  to 
adjust  first  derivatives  of  a set  so  as  to  match  the  subsequent  mean  Kepler 
element  set  at  its  epoch,  and  to  generate  a trajectory  at  any  time  based  on 
the  most  recent  corrected  element  set.  The  original  single  element  set  mode 
is  automatically  effective  when  no  other  elements  follow  since  the  last  element 
set  is  not  adjusted.  Long  term  predictions  beyond  available  sets  are  unreliable. 


Higher  order  derivatives  are  not  adjusted  because  of  their  smaller  effect 
over  a few  days,  and  because  of  uncertainties  over  extended  periods.  For 
significant-drag  orbits  the  rate  of  change  of  mean  motion  is  adjusted  to 
match  the  mean  anomaly;  in  low  drag  cases  the  mean  motion  is  adjusted  instead. 
The  original  and  revised  first  derivatives  can  be  inspected  for  abrupt 
changes,  which  usually  suggests  a transcription  error  in  a data  set.  Plots 
of  satellite  apogee-perigee  such  as  Figure  1 in  Section  5 can  also  point 
up  poor  element  sets  when  the  curves  are  abnormal. 

High  quality  ephemerides  were  prepared  for  the  above  and  other  satellites 
using  daily  elements  covering  periods  of  a month  or  more.  For  reasons 
of  consistency,  earlier  availability,  and  fewer  discrepancies  (see 
Section  4.4),  the  SCF  position-velocity  element  sets  have  recently  been 
favored. 

4.3  Input  of  SCF  Position-Velocity  Vectors 

Program  LOKANGL  has  been  altered  to  predict  satellite  ephemeris  and  look 
angles  given  initial  values  of  position  and  velocity.  Period  decay  (sec/ 
rev)  is  also  required.  Subroutine  DELEM  has  been  added  to  handle  the  re- 
quired conversion  to  SPADATS  parameters 


4.3.1  Computation  of  Initial  Mean  Elements. 

Initial  mean  Keplerian  Elements  are  computed  from  the  input  cartesian 

position  aic  velocities  by  subroutines  taken  from  program  CADNIf^.  Initial 

position  and  velocity  are  assumed  to  be  osculating.  Mean  motion  is  calculated 
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in  revs/day  according  to  Escobal , ' ' 
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i (1-Msin2  i) 
a^l-e^r5 


X 86400. 


where /i  = Earth  gravitation  constant  = 398603.2  — 

sei 

J2  = Second  Harmonic  = .0010827549 
a = Semimajor  axis 
i = inclination 
e = eccentricity 


(2) 
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4.3.2  Derivatives 
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First  and  required  higher  order  derivatives  are  computed  from  secular 
changes  in  the  elements  over  the  first  3 successive  revolutions  due  to 
the  second  harmonic  and  drag.  The  derivatives  of  a function  f are  computed 
as  (M..  is  in  revs/day): 
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where  i£is  the  value  of  f at  the  start  of  the^  revolution. 

4.3.2. 1 Changes  per  Revolution 
a)  Drag 

The  secular  changes  per  revolution  in  the  elements  due  to  atmospheric  drag 

3 

are  computed  according  to  King-Hele 


da  = -a2  9 s (1  + e cos  E)  ^ 
de  (1  - e cos  E)  ^ 


de  = a p 6 


e cos  ^2  (1-e2)  cos  E 


dE 


(1  + e cos  e\ 
1 - e cos  E/ 


(4) 


(5) 


where  E is  eccentric  anomaly,  f is  the  atmospheric  density  and  g is  a 
drag  factor  determined  from  the  input  period  decay.  Changes  in  inclination, 
ascending  node,  and  perigee  argument  due  to  drag  are  currently  neglected 
since  these  are  proportional  to  the  Earth's  rotational  velocity.  However, 
coding  has  been  provided  to  include  expressions  for  these  (lines  19-21 
in  subroutine  INTGND)  if  warranted. 
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The  above  expressions  are  integrated  over  one  rev  (E=0  to  E=2frT>y  Simpson's 
(4) 

rule'  ' using  1 degree  integration  steps.  The  elements  are  assumed  constant 
during  the  integration  of  one  revolution  at  the  values  assigned  to  them 
at  the  beginning  of  the  revolution.  The  density  is  as  given  by  the  DENSEL 
model  using  subroutine  DENSEL  of  CADNIP^. 

b)  Determination  of  6 

The  drag  constant  6 is  determined  by  requiring  the  period  decay  for  the  first 
revolution  to  be  the  input  value,  i.e. 


6 = P/(dP/d6) 


dP  _/86400 

® Ih^ 


(6) 


aa  + 


*5- 


Ai 


where  Aa,  Ae,  and  a 4. are  the  changes  in  the  elements  in  the  first  rev  com 

puted  from/  =1  and  the  partials  of  M with  respect  to  the  elements  are 
computed  directly  from  E^  (1). 

c)  Second  Harmonic 

The  per  rev  changes  due  to  the  second  harmonic  J2  are  as  in  Escobalf2^ 
Ascending  node 

An  = C cos /.  (7) 

Argument  of  Perigee 

A(0  = -c  (2-2.5  sin \ ) (8) 

where 


4.4  Problems  with  ADC  elements. 


4.4.1  Mean  Semi -Major  Axis 

A technical  document^ provided  at  time  of  conversion  from  the  5-card  to  the 
2-card  set  indicated  how  the  mean  (Kozai)  semi-major  axis  a could  be  recovered 
from  the  mean  motion  n and  the  other  data,  since  it  was  now  unavailable. 


where 


r,-  <-3/2i  jAV1  ,1-  ) (i- 


and  <**  = (-3/2)  J2fe2/a,2)  (1‘e  ) [l-(3/2)  sin2i]  (ll) 

>*  = 398603.2,  J2  = .0010827549,  a = 6378.165. 

This  equation  is  in  fact  consistent  with  a and  the  mean  motion  n provided 
on  the  5-card  sets.  Further,  Jan.  1974  2-  and  5-  card  element  sets  for 
satellites  1085  and  1620  are  effectively  identical,  indicating  no  disparity, 
due  to  the  transition. 

Investigation  of  first  order  secular  perturbation  theory^  however  discloses 
that  the  Kozai  mean  semi-major  axis  is  improper  for  use  as  the  mean  semi- 
major axis  in  programs  L0KANGL  and  R0PP.  The  oblateness  of  the  earth  is 
considered  to  perturb  the  mean  motion  but  not  the  semi -major  axis  a. 

If  unperturbed  n,  = ./'/' 3 
Then  perturbed  n = n„  (1-gC  ) 

As  before,  defining  a*=  % ) ^3 

We  have  a = a0  Jmi/3)«£  - (l/9)£2j  .(13) 

Qualitatively,  for  inclinations  near  90*  ^is  positive  and  the  perturbed 
mean  motion  n is  lower  than  n#  . Hence  the  fictitious  atf  must  be  larger 
than  the  mean  semi -major  axis  a. 

Figure  3 shows  before  and  after  printer-plots  of  the  relative  mean  semi- 
major  axes  as  obtained  from  ADC  (A)  and  SCF  (X)  elements.  Use  of  Equation 
13  instead  of  Equation  9 makes  the  values  from  the  two  sources  consistent. 
4.4.2  Irregularities 

Figure  4 shows  inclination  obtained  from  ADC  and  SCF  elements  over  a 60 
day  span.  Results  of  R0PP  prediction  runs  using  C^A  adjusted  to  give  a 

HT 

matched  initial  nodal  period  are  also  included  . A similar  plot  for  ec- 
centricity was  shown  in  Figure  4 of  Section  3.  It  can  be  seen  that  the  ADC 
mean  elements  are  much  more  erratic  than  the  SCF  derived  values.  The  sys- 
tematic sinosoidal  difference  in  eccentricity  seen  in  Figure  3 of  Section  3 
occurs  because  the  ADC  elements  provide  mean  mean  eccentricity,  and  the 
adjustment  for  the  J3  zonal  is  needed  according  to  Equation  1. 
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Figure  4.  ADC,  SCF  and  ROPP  (predicted)  Inclinations 


Figure  5.  Comparison  of  ADC  and  SCF  Eccentricities  After  J,  Correction  to  ADC  Elements 
See  Section  3,  Figure  3 for  pre-correction  plot. 
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Figure  6.  ADC,  SCF  and  ROPP  (predicted)  Perigee  Altitudes 


Figure  5 compares  the  ADC  and  SCF  eccentricities  after  the  above  correction 
has  been  applied.  A small  systematic  difference  is  still  present,  and  the 
irregularities  are  evident.  Figure  6 compares  the  corresponding  perigee 
altitudes  over  a 4 month  period.  The  SCF  values  are  seen  to  be  consistently 
smoother.  Incidentally,  the  lower  density  scale  height  is  evident  in  ROPP 
from  the  faster  altitude  decay. 

4.5  Other  Modifications 

Program  LOKANGL  was  compacted  by  elimination  of  unused  routines,  variables, 
and  some  double  precision  code.  A number  of  features  were  also  added  and 
are  described  below.  Card  input  format  has  been  modified  minimally  and  is 
shown  in  Figure  7. 

4.5.1  Output  of  Osculating  Position-Velocity. 

Routines  taken  from  CADNIP^  convert  mean  Keplerian  elements  to  osculating 
elements  which  include  the  short  period  perturbations.  Position  and  velocity 

are  then  directly  obtained.  Keplerian  elements  when  requested  are  output  as 
mean  values. 

4.5.2  Binary  Ephemeris  File 

The  binary  output  file  (TAPE3)  is  identical  in  format  to  the  file  generated 
by  DABOS^'  except  that  the  maximum  number  of  stations  has  been  limited  to 
20,  and  the  observational  data  for  all  stations  (up  to  20)  are  included 
with  the  ephemeris  record. 


4.5.3  Sub-ionospheric  Option  in  LOKANGL. 

Subroutine  SILL2  was  optionally  auoed  to  program  LOKANGL  in  order  to  cal- 
culate azimuth,  elevation,  longitude,  and  geodetic  latitude  of  the  location 
below  a given  ionospheric  height  along  the  direction  to  the  satellite. 
Subroutine  SILL2  was  provided  by  the  Ionospheric  Physics  Laboratory. 

4.5.4  Program  LOKANGL  has  undergone  a number  of  changes  in  order  to  adjust 
and  present  the  input/output  information  in  compact  form.  Also  the  subrou- 
tines WRTAP  and  WRSTA  which  generate  observational  printout  were  unified 
into  one  subroutine  WRSTP  which  allows  optional  output  by  time  or  by  station 
only. 


— 
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INPUT  FORMAT  FOR  PROGRAM  LOKANGL 


OATA  DECK  SETUP 
CARO  COLS 


DESCRIPTION 


(jwVJL  wT 

tmu 

IPAfc 

RHINA  t tk 

W t-OKft 

l=FORH 

NO. 

3 \ 

SCF 

2-CARO 

POS.-VEL.  SET 

2- FORM 

NO. 

2 \ 

ADC 

2-CARO 

ELEMENT  OATA  SET 

5*F0RH 

NO. 

1 \ 

AOC 

5-CARO 

OATA  SET 

ONE  OR 

MORE 

ELEMENT 

SETS 

1-3  NO.  OF  STATIONS.  IF  NO  STATIONS  USE  0 
6 COOE  1 OR  0 FOR  PRINT  CONTROL  OF  STATION  DATA 
••PRINT  0Y  TINE  ONLY 
1*PRINT  BY  STATIONS 


I ONOHT*SUB- IONOSPHERIC  ALTITUOE  (KM)  13 

STATION  LOCATION  CAROS. IF  NMS.GT.O 
1-5  I.O.  OF  STATION  (NUMBtRi 

8 COOE  2=G£00ETIC  SYSTEM 

0- 23  STATION  G£06€TIC  tATITUOE.  DEGREES 

24-38  STATION  L0N6ITU0E  (POSITIVE  WEST).  DEGREES 
39-53  STATION  HEIGHT.  METERS 
61-72  NAME  OF  STATION 

1- 15  TIME  INCREMENT  IN  SECONDS.  F15.5 


START  COLS\  17-18? 20-21 ? 23-24? 26-29 ? 31-34? 36-41 

PINAL.  CGLSV  4-3— 94?  46- 47?  49- 50-?  52-55  ? 57-69 ? 62  -67 

TIME  MONTH?  DAY?  YEAR?  HOUR?  MIN?  SEC 

FORHAT  12?-  12?  12?  F*.l?  F*.l?  F$.3 


2 IMPOSITION  AND  VELOCITY 

-4 1=SU6- SATELLITE  DATA 

6 1=ME AN  KEPLERIAN  ELEMENTS 

-8  t=09S€RV  ATICN  OATA 

9-10  TAPE  3 CONTROL?  STANOARO=0 

11-13  STANOARO-l.O 

1-14  "END  OF  PROBLEM- 


ABOVE  MAY  BE  REPEA TEO  N TIMES 
I.AST  CAROX  J!9S!-IN-IJOL.  -t 


Figure  7.  Card  Input  Description  for  Program  LOKANGL 
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5.0  Satellite  Orbit  Support  Programs 

Initiator:  I.  M.  Hussey 

Project  No.:  0001  Problem  No.:4506,  4801 


5.1  Apogee- Peri gee  Program  AP 

In  order  to  graphically  present  and  evaluate  various  satellite  orbi ts,  program 
L0KANGL  is  often  executed  simultaneously  with  program  AP  which  calculates  apogee/ 
perigee  data  using  the  minute-by-minute  binary  ephemeris  file  TAPE3  generated  by 
program  L0KANGL.  Program  AP  has  undergone  changes  which  allow  simultaneous  cal- 
culations of  apogee  and  perigee  data  which  permits  more  efficient  use  of  the  system 
of  programs  L0KANGL,  AP  & PL0TIT  (Section  5.2).  Poor  quality  element  sets  or 
transcription  errors  are  readily  identifiable  as  irregularities  in  the  apogee- 
perigee  plots. 


Initiator:  I.  M.  Hussey 

Project  No.:  0001  Problem  No.:4506,  4801 

5.2  Program  PLOT IT 

Program  PL0TIT  performs  three  different  types  of  plotting.  They  are,  1)  indi- 
vidual, 2)  summary,  3)  multiple  plots. 

Individual  plots  have  one  curve  only  per  graph.  Summary  plots  contain  several 
curves  of  the  same  variable,  ear*  cu^ve  representing  the  use  of  a different  data 
set.  Multiple  plots  contain  several  curves  of  different  variables  from  the  same 
data  set.  The  program  was  originally  developed  by  C.  Foley  (Boston  College) 
and  has  undergone  a number  of  changes,  a)  Option  to  read  variable  data  from  any 
tape  (not  only  from  punched  cards)  was  added,  which  allows  execution  of  PL0TIT 
in  conjunction  with  other  programs  e.g.  L0KANGL,  AP.  b)  Appropriate  subroutines 
of  PL0TIT  were  modified  in  order  to  include  satellite  revolution  number  on  the 
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*****  OArA  DECK  SET J* 


NPASS 


THRJST  CARDS  A RE  ENTERED  NEXT  AS  FOLLOWS... 

M "t  1 tutT^  v/r  f *Ttt7  j t uvCuK 

ENCE 


4-3 

-11-53- 

51-52 

p i ^ ^ 


53-61 


66-71 

73-74 


-1-mOICATES  END  OF  THR.)  12 


VARIAB.E  ABBREVIATION  A5 

-4XIS  L-A5EI  -OR- VARIABLE MH 

VAR.  STM.  CODE  NlM  BE R 12 


OWER  -IMIT  OF  VARIABL 


T AXIS  HGT.  IN  INI  HES  F4.0 

NUMBER  OF  DIGITS  TO  RT.  OF  .12 


OLS  54-56  ARE  FIX,  THEN 


3ER  OF  DIGITS. 


****************************************** * ** * *»*,«*** ********< 


AS  MANT  AS  15  MORE  mH>Zm  CAROS,  ONE  3ER  VARIABLE. 

■»  r»  T u ki  . naa  ^ A nn  f*  r .*«»  ..  • a 


THEN  THE  OAT  NUMBER  CX-AXIS)  IS  CA.CULATED  ASSUMING  THAT  THE 


NUMBER,  2.)  THE  HOUR,  AND  3.)  THE  MINUTE. 


ARDS) 


MULTI? (I, K)  CONTAINS  THE  VARIABLE  FIELD  NUMBER.  IT  ALLOWS  UP  TO 

6 MJ-TI^LE  PL3TS  EACH  HAVING  AS  HAN  AS--5  VARIABLES— (RT.  JUSTI 

FIED) 


12  -IELD)  WIL.  INSURE  THE  “LOTTING  DF  THE  VARIABLE  WHOSE  FIELO 


• 

N+4 

1-32 

i -An 

NCASES, ISUMRV(I) . 

Mill  TTO  / V t/t 

1612 
tn  verii 

* 

• 

• 

N«-5 

A 3U 

1-32 

"UL  f X“t  X f" 

-15~,I?LT(I) 

O !5iC> 
1612 

69-70  VARIAB.E  DATA  FIL:  NUMBER  (1=CARD  INPUT) 


ONES  THE  VARIABLE  DATA 

— * ft?  4A-T-FI 3 ,-312,  *-2ifl  , 2F5.  1 , F 7. 2,  ?£9rTjF7 .2,E9»3,F6§3, 

F4. 1, I2»F4. 0 , 14) 

-T-HI 5— 3-ATA— MUST— 8E-F-QLL-DWF0— B¥— AN-EH-D— OF— FILE, — AND— TH€- 
SATEL.ITE  AND  MODEL  USED  CARD,  AND  VARIABLE  DATA 


Figure  1.  Card  Input  Description  for  Program  PLOTIT 
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horizontal  scale,  as  well  as  to  allow  correct  numbering  of  the  horizontal  axis 
for  day  number  over  consecutive  years,  c)  Thrust  times  can  be  indicated  on  the 
plots  (see  Figures  in  Section  7). 

Program  PLOTIT  was  regularly  used  in  combination  with  a number  of  different  programs 
(ROPP,  LOKANG,  AP)  Latitude  and  local  time  plots  for  apogee  and  perigee  of  various 
satellites  were  generated. 

For  reference,  the  data  deck  setup  is  described  in  Figure  1.  Figure  2 shows 
separate  perigee  and  apogee  plots  that  were  generated  for  satellite  8468  (S3-2) 
using  program  ROPP. 


Initiator:  I.  M.  Hussey 

Project  No.:  0001  Problem  No.:  4801 

5.3  Satellite  Visibility  Program  MINVIS 

Given  time  span,  satellite  elements  and  epoch,  and  observer  coordinates,  program 
MINVIS  calculates  minutes  of  optical  visibility  per  day  at  each  station. ^ 

Output  is  in  tabular  form  and  also  as  a histogram  per  station.  Visibility  also 
requires  the  sun  angle  to  be  within  specified  limits. 

This  program  was  converted  for  use  on  the  CDC6600  under  Fortran  Extended.  Minor 
corrections  were  made  and  the  program  was  checked  out.  Card  input  format  is  as 
follows: 


(1) 


Kellaher,  J.G.,  "Satellite  Observability",  Analysis  and  Simulation  Branch, 
Problem  1444. 


Card  Input  Format  for  Program  MINVIS 


CARD 

VARIABLE 

FORMAT 

1 

Run  No. 

X3 

Begin  YR,  MO,  DY,  HR,  MN,  SC 

612 

End  YR,  MO,  DY,  HR,  MN,  SC 

612 

At  (sec) 

15 

2 

Satellite  Ident. 

3A6 

3 

Epoch  of  Elements 

YR,  MO,  DY,  HR,  MN,  SC 

612 

4 

Elements 

Semi-axis  (km) 
Eccentricity 
Inclination  (deg) 

Mean  Anomaly  (deg) 
Arg.  of  Perigee  (deg) 
Asc.  Node  (deg) 

6E13.8 

5+ 

Station  Cards  (up  to  60) 

Station  Ident. 

A 6 

Station  Network 

11 

Station  No. 

Z3 

Latitude  vd=g.,  min.,  sec.) 

5X.3F5.0 

E.Long.  (deg.,  min.,  sec.) 

3F5.0 

Altitude  (meters) 

F7.0 

Minimum  Sun  Elevation  (deg) 

F4.0 

Maximum  Sun  Elevation  (deg) 

10X.F6.0 

Last 

"END" 

A3 
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6.0  Satellite  Orbit  Determination  Prog ram- DABOS 


Initiator:  I.  M.  Hussey 

Project  No.:  0001  Problem  No.:  4801 

(1  2) 

Program  DABOS  v ’ ' is  unique  compared  to  other  programs  used  at  AFGL 

such  as  CELEST,  TRACE  and  CADNIP  which  conventionally  optimize  initial 

conditions  to  fit  extensive  models  of  the  earth's  geopotential  and 

atmosphere  to  radar  tracking  or  Doppler  signals.  In  DABOS  a Kalman 

filtering  procedure  is  used  which  dynamically  weights  observations  vs. 

model  predictions  along  the  trajectory  in  both  forward  and  backward 

integration  modes.  This  procedure  is  particularly  desirable  when  good 

quality  observations  are  available  and  the  model  is  relatively  poor. 

On  the  other  hand  a sequence  of  marginal  observations  from  one  station 

can  cause  offsets  to  build  up  in  the  trajectory,  leading  to  erratic 

apogee-perigee  altitudes  and  improper  rejection  of  data  elsewhere  in  the 

trajectory.  A hybrid  Kalman  filtering  and  smoothing  approach  had  been 

initiated  but  not  debugged  by  the  original  investigator.  Refinements 

were  investigated  including  dynamic  drag  adjustment  on  low  perigee 

satellites,  smoothing  over  extended  sets  of  observations,  and  integration 
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of  the  station  bias  estimation  program  BIAS'  ' into  DABOS.  The  latter 
was  converted  for  operation  on  the  CDC  6600.  Other  tasks  that  were 
completed  allow  input  of  SCF  observational  data  in  "DECOR"  format,  and 
incorporation  of  advanced  geopotential  and  atmospheric  models  in  DABOS. 

Another  valuable  aspect  of  DABOS  is  its  output  and  ephemeris  tape 
compatibility  with  program  L0KANGL  which  itself  was  originally  derived 
from  DABOS.  This  plus  its  relatively  modest  operational  requirements 
(110K  with  two  overlays)  permit  flexibility  in  development  and  research 
investigations.  Parallel  aeve < client  of  these  programs  has  promoted 
retention  of  these  attributes. 

6.1  SnflSothinq 

Some  difficulty  arises  in  the  use  of  program  DABOS ^ over  large  time 
spans  in  which  data  is  lacking  and/or  high  geomagnetic  activity  occurs. 

For  example,  large  adjustment  in  the  state  vector  by  the  Kalman  filter 
for  the  first  few  data  points  after  such  time  spans  occasionally  produces 
noticeable  jumps  in  perigee  heights  (l-2km).  A recent  version  of  DABOS 
which  incorporates  an  ad-hoc  smoothing  technique^was  therefore  tested 
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and  found  useful  in  alleviating  these  discontinuities.  Changes  and 
usage  of  this  program  are  described  here. 

6.1.1  Smoothing  Program 


The  smoothing  version  of  DABOS  uses  the  filtered  state  vectors  at  6 
data  points  to  refine  the  vector  at  the  4th  as, 

u (t4>  -I  2‘  lk'4*  W/1 2'  U) 

k=l  k=l 


where  X^t^)  is  the  state  vector  at  t^  obtained  by  integrating  from  the 
filtered  state  vector  at  t^  to  time  t^.  Consecutive  state  vectors  in 
Eq  ( 1 ) must  be  separated  by  at  least  10  minutes  in  time.  The  Keplerian 
elements  rather  than  position-velocity  are  used  in  this  procedure.  The 
resulting  smoothed  vector  is  used  for  subsequent  orbit  prediction  (in  the 
filtering)  but  the  original  filtered  vectors  (not  the  smoothed  vectors) 
are  saved  for  subsequent  smoothing. 

An  additional  feature  of  the  smoothing  version  of  DABOS  is  the  calcu- 
lation of  the  state  transition  matrix  by  the  "perturb  and  integrate" 
method: 


, xk  (X,  ^X,)  ■ Xt  (X.  ) 

ax,  ax, 


where  X -j  is  a component  of  the  state  vector  at  time  ti  , and  the  state 
vector  X^(X^  ) is  obtained  by  numerical  integration  from  ti  to  t^. 

6.1.2  Changes 

For  addition  of  atmospheric  density  subroutine  DENS  see  Section  6.3. 

For  expanded  geopotential  models  see  Section  6.4. 

For  addition  of  DECOR  input  and  other  changes  see  Section  6.5. 

6. 1.2.1  Restricted  BCD  ephemeris  output 

A new  input  parameter  NPRT  has  been  added,  specifying  the  ratio  of  BCD  to  binaiy 
ephemeris  output  interval  when  both  BCD  and  binary  output  are  requested. 

This  permits  a fine  interval  for  the  binary  ephemeris  when  such  is  necessary 
for  further  processing  such  as  for  perigee  computation,  while  at  the  same 
time  avoiding  an  overly  cumbersome  BCD  printout. 
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6. 1.2.2  Variable  number  of  points  in  smoothing 

Coding  has  been  added  to  vary  the  number  of  data  points  used  in  smoothing. 
In  addition,  smoothing  can  be  suppressed. 

6. 1.2. 3 Atmospheric  Density  Update 

The  smoothing  version  of  DABOS  was  tested  using  2,4,  and  6 smoothing 
points  and  found  to  be  reasonably  successful  in  removing  perigee  height 
jumps  found  previously.  This  however  does  not  establish  the  accuracy 
of  the  resulting  orbits.  It  is  evident  that  the  real  cause  of  the  orbital 
discontinuities  encountered  here  are  dynamic  variations  not  properly 
accounted  for.  The  smoothing  technique  will  average  over  these  variations 
but  not  reproduce  their  details.  Furthermore,  a DABOS  - type  run  which 
includes  atmospheric  density  updates,  along  with  state  vectors  could 
produce  a record  of  time  variation  of  density  which  could  serve  as  a 
useful  research  tool.  For  these  reasons  it  was  decided  to  introduce  a 
coefficient  C which  multiplies  the  model  density. 

The  procedure  for  updating  this  parameter  is  similar  to  that  given  on 
pages  32-4  of  Ref.l  for  the  determination  of  type  1 (dynamical)  con- 
stants. However  modification  is  necessary  to  account  for  the  fact  that  the 
parameter  to  be  determined  here  is  not  constant.  This  is  done,  in  the 
calculation  of  the  predictive  covariance  matrix,  by  adding  to  the  density 
covariance  a term  Pc"  which  represents  the  uncertainty  in  density  varia- 
tion since  the  previous  observation.  Thus  Eqs.  (60)  and  (62)  are  modified 
to  read 


')  - D*  + p ft 

c C C 


P*  is  currently  treated  as  input  with  a value  of  .05  found  to  be  most 
satisfactory. 

6.2  Observation  System  Bias  Estimation 

E»o*r1ence  with  satellite  tracking  has  shown  that  observations  usually 
mam  balses  or  systematic  errors  in  addition  to  the  random  noise. 

« tracking  systems  (particularly  electronic  systems),  two  more 

»<•  m xnplicate  the  situation.  The  ratio  of  the  random  noise  to  the 
*«  usually  small,  and  the  functional  (often  time  dependent) 
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form  of  the  bias  errors  is  almost  impossible  to  predict  and  formulate. 

This  Is  not  due  to  any  deficiencies  in  theory  but  simply  to  the  fact  that 
equipment  can  malfunction  and  often  does  in  a manner  that  cannot  be  pre- 
dicted or  even  recognized  from  the  output.  On  the  other  hand,  high  fre- 
quency random  noise  and  sudden  large  changes  in  the  output  parameters 
can  be  easily  recognized  and  treated  accordingly. 

In  the  orbit  prediction  program  DABOS,  this  problem  was  alleviated  by 
estimating  the  system  standard  deviations  inclusive  of  the  bias  errors 
by  successive  approximations  using  the  minimum  variance  method  in  a 
multiple  filtering  process.  While  this  method  has  certain  advantages  and 
might  be  the  only  practical  method  in  many  cases,  it  might  not  give 
sufficiently  accurate  estimates  of  the  orbit  in  case  the  bias  errors 
are  very  large  and  the  number  of  observing  stations  is  small. 

It  is  possible  to  go  one  step  further  and  try  to  estimate  the  biases 

simultaneously  with  the  orbital  elements.  A development  based  on  the 

minimum  variance  method  and  outlined  in  Ref.  1,  pp  34  to  37  is  utilized 
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in  the  present  program  BIAS'  ' . This  program  was  reviewed  and  debugged 
for  operation  on  the  CDC  6600. 

The  program  was  found  to  yield  satellite  positions  deviating  from  those 
of  the  original  DABOS  by  up  to  50  km,  and  station  location  biases  of 
the  same  magnitude.  These  results  have  been  attributed  to  incorrect 
calculation  of  the  bias-orbit  covariance,  and  this  was  corrected. 

6.3  Atmospheric  Density  Models 

An  objective  of  acquiring  and  operating  program  DABOS  was  to  facilitate 

atmospheric  density  studies  through  the  calculation  of  the  atmospheric 

drag  on  satellites.  To  further  this  end,  several  density  models  have 

been  added  to  DABOS  as  optional  alternatives  (including  the  Jacchia  1964). 

This  was  done  by  adding  subroutine  DENS  which  had  been  written  for  program 
(5) 

CADNIP'  , another  orbit  determination  program.  Subroutine  DENS  contains 
the  option  to  compute  the  density  according  to  any  of  the  models  listed 
in  the  Input  Format  description  at  the  end  of  this  section.  (The  missing 
numbers  correspond  to  other  density  models  handled  by  CADNIP,  but  not 
in  this  subroutine.)  These  models  contain  explicit  dependence  on  many 
of  the  following  effects: 

(a)  Solar  activity 

(b)  Geomagnetic  activity 

(c)  Diurnal  effect 
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(d)  Semiannual  effect 

(e)  Polar  bulge 

(f)  Seasonal -latitudinal  effects. 

The  density  is  computed  by 

log  1Qf  = log10  f(T^  h)  + g 

where  h is  the  satellite  height,  log^f  is  computed  from  a table  supplied 
as  input,  and  g and  T*  (exospheric  temperature)  are  computed  analytically 
as  functions  of  the  above-mentioned  effects. 

Upon  an  initialization  entry,  the  requested  density  model  table  is  read 
and  the  appropriate  section  of  the  solar  and  geophysical  file  is  read. 
Upon  a regular  entry,  the  input  position  is  converted  to  right  ascension, 
declination,  and  height.  The  appropriate  formulas  for  computing  local 
exospheric  temperature  are  exercised.  Two  dimensional  quadratic 
interpolation  in  temperature  and  height  yields  the  desired  density  which 
may  be  corrected  for  geomagnetic  index,  seasonal  latitudinal  variations, 
or  semi-annual  variations. 

Program  DABOS  currently  inputs  the  following  values  to  DENS: 
Initialization  Inputs 

EQTRD  -equatorial  radius  (6378.165) 

PLRRD  -polar  radius  (6356.783) 

TGRMJD  -reference  time  in  M.J.D.  (40222.0) 

THET7R  -sidereal  angle  of  Greenwich  at  reference  time  in  revolutions 
(.27907816) 

THETDT  -modified  rotation  rate  of  Greenwich  in  revs/day  ( .27379093E-2) 
MODTYP  -desired  model  number 

TAPE  8=NDENTP  -logical  number  of  file  containing  atmospheric  table 

TAPE  9=NAPFTP  -logice.1  r>'..r,her  of  file  containing  solar  and  geomagnetic  data 

TSTAPF  -time  in  M.J.D.  to  at-rt  reading  and  accumulating  solar  and 
geophysical  data  = DJULN+33282 

Regular  Inputs 

TIME  - time  in  M.J.D.  of  desired  density  = DT0TSP+33282 
PV  - position  in  normalized  xyz  inertial  cartesian  coordinates 
DENFAC  - multiplicative  factor  to  apply  to  densities 
Program  DENS  outputs  the  following  values  to  subroutine  MOTION: 

Outputs 


TE  - exospheric  temperature 


HEIGHT  - height 

DENSIT  - density  in  cgs  (g/cc) 

The  significant  part  of  new  density  facilities  for  program  DABOS  is 
subroutine  APFTMP  which  computes  array  APFTAB  of  19  x 200  values  for 
storing  200  days  of  solar  and  geomagnetic  activity  data  and  also  array 
FBAR  of  308  values  of  F10.  7.  This  subroutine  has  undergone  changes 
which  allow  calculation  of  the  above  tables  in  case  of  early  termination 
of  the  solar  flux  file.  In  this  case  unless  the  "insufficient  data" 
message  is  printed  the  insufficient  data  is  replaced  in  averaging  cal- 
culations by  the  last  record. 


6.4  Earth's  Gravitational  Field.  (GF0RCE  addition) 

The  mathematical  model  of  the  dynamical  system  is  expressed  by  the  equa- 
tions of  motion  with  terms  representing  the  various  forces  acting  on  the 
satellite  including  Earth's  gravitational  field,  atmospheric  drag.  Sun's 
and  Moon's  gravitational  fields  and  solar  radiation  pressure. 

6.4.1  Gravitational  Model 

The  mathematical  representation  of  the  Earth's  gravitational  field  is 
expressed  by  means  of  gravitational  potential  function  which  can  be 
written  as 


(C  cos  mA  + sir 
v nm  nm 


C = -J  = C , S =0 
no  n n - no 

where^/C  = G Me 

J , C . and  are  numerical  coefficients,  Rr  the  mean  equatorial  radius 
n nm  nm  E ^ 

of  the  Earth,  r the  distance  of  the  satellite  from  the  center  of  the  Earth, 
0the  latitude,  and  Pn  the  associated  Legendre  polynomial 

m,x)  = d.^m/S  (X) 


where  Pn  is  the  Legendre  polynomial.  The  longitude^is  to  be  counted 
positive  to  the  east  in  this  application. 

The  harmonics  represented  in  the  gravitational  potential  function  are  called 
spherical  harmonics.  If  0<m  < n they  are  called  tesseral  harmonics  as  a 
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special  case  of  the  spherical  harmonics.  If  m = 0,  they  are  called  zonal 
harmonics,  and  if  m = n,  they  are  called  sectorial  harmonics.  The  gravi- 
tational potential  for  bodies  with  spherical  symmetry  can  be  expressed 
by  the  zonal  harmonics  only,  i.e.,  the  potential  is  a function  of  latitude 
and  independent  of  longitude.  For  bodies  of  arbitrary  shape,  the  potential 
must  include  the  tesseral  harmonics,  which  are  dependent  on  both  latitude 
and  longitude. 

6.4.2  Addition  of  GFORCE  to  DABOS 

In  order  to  examine  various  geopotential  models  and  to  extend  the  capa- 
bilities of  DABOS  subroutine  GFORCE  which  includes  calculation  of  tesseral 
harmonics  was  adapted  from  program  CELEST.  The  SEIII  and  GEM  models 
were  tested  against  the  single  original  model  in  DABOS,  which  included 
calculation  of  only  the  first  six  zonal  harmonics  of  gravitational  potential. 
For  these  large  g.p.  models  the  execution  time  increased  by  well  over  an 
order  of  magnitude. 

Geopotential  models  for  use  by  DABOS  should  be  created  on  TAPE25.  Format 
is  .similar  to  that  used  for  input  to  CELEST: 


Variable 

Format 

Card  1 

Alphameric  I dent 

No.  of  Cards  following(N) 

p 

mmax  ’ max 

A10 

£20.14 

-2E20.14 

Cards  2 

m,  JL 

5X,2I5 

to  N + 1 

c s 

term’  term 

2E20.14 

Note:  If  N is  negative,  print-out  of  the  geopotential  model  terms  is 
suppressed. 


6.5  Other  Changes 

Program  DABOS  was  changed  in  order  to  allow  input  of  SCF  observation  data  cards 
in  "DECOR"  format.  Option  to  read  cards  in  this  format  was  also  added  to  program 
BIAS  and  the  smoothing  version  of  DABOS. 
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Program  DABOS  was  also  modified  to  allow  input  of  pre-launch  elements  in  the 
same  format  as  the  Keplerian  or  the  position-velocity  elements. 


The  following  pages  describe  the  Input  Format  in  effect  following  these  various 
changes. 


p 
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INPUT  FORMAT  FOR  PROGRAM  DA30S 


CO  1 


TITLE 

-ANY  DESCRIPTIVE  TITLE 


HI  S C .INFORMATION 

1-5  ND  SAT  SATELLITE  NO. 


i=JACCHIA  1964 
2=1966l  SUPPLEHENTS. 

3= JACCHI A 1971 
4=1968  CHAMPION 
5=JACCHIA-WAL<ER-BRUCE 


rnsn 


c 

_c 

COL  11 


7=LOCKHEEO/NASA 

&=  CH  ANDRA- KRISHNA  MURTH  t 

12=J  ACCHI A 1 970 


1=REJECT  OBS.  WITH  ELEVATION  NOT  WITHIN 
LIMITS  SET  ON  CARO 


0=00  NOT  REJECT  ANY  OBSERVATION  CARO 


1= CORRECT  ELEVA' 

0=NO  CORRECTION 


RATE.  FOR  REFRACTION 
13 


=N0  INPUT 


COL  26 
C 


G=N0  PRINT 

26  ISMOOH  1=APPLY  SMOOTHING 

M 


28-29  NSTPRT  NO.  OF  STATIONS  IN  EPHEMERIS  PRINTOUT 


jJiiT  E_ 


(IK 


2.0=OOUBLE  CORRECTION 


COL  63-65  IlUi  REVOLUTION  NO.  AT  TIME  OF  INPUT  ORBIT AL  ELEMENTS 


COL  66-74  TIMTO.  SYSTEM  TIMING  ERROR  (SEC.) 

COL  76 I 0_R ELEMENT  TYPE  INPUT 

C 0 =EL EMENTS  INPUT 


F9.3 


iU'KI 


OKI 


i'lUUI 


C 

JC 

CO  3* 


-1=PRE  LAUNCH  INPUT 


CO  3+  STATION  OATA  CAROS  ( NMS  CAROS) 

GQ.L  E-.5  . NUilSTA  STATION  NUMBER  

COL  10  N3G  1=EL.  A AZ.  IN  GEOCENTRIC  SYSTEM 


COL  16-31  PHILAT  STATION  GEOOETIC  LATITUDE  (DEGREES) 

COL  . 31-45  S- ON  STATION  WEST  LONGITUDE  (DEGREES) 

COL  46-60  A.T  STATION  ALTITUDE  (M.) 


COL  4 

U 

CO  3+  + 


E15.8 

E15.8 

E15.8 


END  OF  STATION  OATA 


C 

CO  4 
COLUMNS 


1-15 


SATELLITE  INFORMATION 

16-31  31-45 


Figure  1.  Card  Input  Description  for  Program  DABOS 
1 of  4 
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COLUMNS  2-3 


C. FORMAT 


F3.1 


CO  6 


ORBITAL  ELEMENTS 


I 0R=-1 


COL  14-2&  ECCEN  ECCENTRISITY 


rJiw  rc  li 


Y3OS  (KM) 


iiJiL€rcii 


r.»  Am  llKtT«l* 


RAOIUS (FT) 


lilf.WnBViafiNi 


COL  40-52  WA  SC  RIGHT  ASCENSION (DEG)  XVELCKM/SEC)  L DNGITUDE( RAD) 


M'iffWiUa'Jirtaaiil 


i l'laBl'CriMaiiW.vjfMhivvri.V  til 


COL  66-75  XMEAN  MEAN  ANOMALY (DEG)  2VEL  (KM/SEC)  FL.  PATH  ANGLE (RAD) 


CD  7 


ftian.u'Ki 


C.  VAR 


INITIAL  ESTIMATE  OF  POSITION  l VELOSITY  ERRORS 


PMAT<1,1) 6E13.7, 


>PMAT(6»6) 


CO  6 


C.  VAR 


MEASUREMENT  ERRORS  (STANDARD  DEVIATIONS) 


QP<1)  3P  ( 2)  Q ( 1) 9-8.4 


,Q(7) 


COLUMNS 


9-16  17-24  25-32  33-49  41-48  49-56  57-64  65-7? 


CO  10 


COLUMNS 


ihikii  a*  11 


C.  VAR 


CO  11 


ELEVATION  ANGLE  ALLOWED 


ihiiiii.  i: 


1-15 


16-30 


IIl  VSIU  l lllil 


lala'Milil' 


EL E MIN 


ELEMAX 


COEFFICIENTS  OF  ZONAL  HARMONICS 


COLUMNS 


16-30 


31-45 


46-60 


61-75 


USE  IF  COLUMN  17  ON  CARO  2 EQUAL  1 


C.  VAR 


E15.8 


COLUMNS 


17-18 


20-21 


23-24 


26-29 


31-34 


36-41 


CONTENT 


MONTH 


YEAR 


HOUR 


MINUTE  SECOND  -START 


C.  VAR 


JM3PRT  JDYPRT  JYRPRT  HRPRT 


XMIPRT  SECPRT  (2) 


Figure  1.  Card  Input  Description  for  Program  DAB0S 
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CONTENT 

COl 6 


COL  73 


MONTH 

M5SPR 


JTYPRT 


YEAR  HOUR  HINUTE  SECOND  -FINISH 

1 1 1 1 a i i 


ujajM 


■Uiki  I ■ ill;  Pi'iaiiiH  i m fli  ■ 1 i fiiMiiikijm  ifiimii  §■ \ 


TYPE  OF  EPHENERIS  PRlNrOUT 
ft  =m n fdmfmpjtc 


1=PRINT  POSITION, VELOCITY 


3=PRINT  POSITION,  VELOCITY, ELEMENTS, STATION  OBS.  OAFA 


5=PRINT  STATION  OBSERVATION  DATA 


-l*BCO  TA«*E  ONLY 


1=BINARY  TAPE  ONLY 


IN9UT  IF  NOSPRI  (COL.  7 
F-e.  a — a 


JSPOA  JSPYR  SPHR  SPHI  SPSE! 


SPSE. 


iT:  VL*  :1T»  i fiTTBiYVm.' 


NERASE 


COL  10-il  NOBS 


COL  15-16  ZO AT 


SATELLITE  NUMBER 


YEAR  LAST  2 DIGITS  OF  19XX 


HOUR  OF  ZULU  TIMES  OF  19XX 


F2.0 


IV4II1 1 


v.irniiwiviiiiii 


COL 

19-23 

ZO  AT 

SECONO  OF  ZULU  TIME 

F5.3 

1 

■hiTI 

?fc-29 

QRSN 

b 

COL 

31-32 

ERASE 

AZIMUTH  OR  RT.  ASC. 

F2.0 

LLVc<:V12 


COL  35-37  ERASE 


COL  U6 


COL  56-63  03SN 


mm 


AZIMUTH  OR  RT.  ASC. 
RANGE 


EXPONENT  (RANGE  =R ANGE*1 E. **IEX) 


ELEVATION  RATE  (OVERPUNCH  0<> 


F3.1 

FT  . F 


II 


F5.L 


COL  68-72  0 3 SM 


1 = EL  E V ATION  AND  AZIMUT 


F5.5 


UxEBImI  H • i 2 1 1 1 • I • I ; IV  * 1 ? II  Hi  r M fl  fi;  I WM  • ici  H3  J Cv  ■•T;l 


COL  76 
r 


CD  15+* 


: vvi*.  rnm’t 


IP  I- ■.HI  H 1 H IA4ii:M:nr.i:] 


COL  15 


NJYR 


3=ELEVATI0N,  AZIHUT,  RANGE,  RANGE  RATE 

l*  = PLF«»TTnN.  A7TMSIT.  PA  Mf.F.  Plur.P  PATF.  I 


RATE,  AZIMUTH  RATE,  RANGE  ACCELERATION 


YEAR  or'  r.  TIT  NO  X 


OBSERVATION  CARDS,  DECOR  FOMAT 

VFHTr.l  F MfMPFP 


ANTENNA  (STATION  TYPE) 


UNITS  OF  YEAR 


NTP=2 


NTP=4 


NTP=7 


mamai 


ON  (DEGREES) 


12 

— FS*J 
F6. 3 


Figure  1.  Card  Input  Description  for  Program  DABOS 
3 of  4 
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7.0  Atmospheric  Density  Studies 


Initiator:  K.  Champion 

Project  No.:  6690 


Problem  No.:  4547 


Program  CELEST  is  a satellite  trajectory  data  processing  program  developed  by 
Naval  Surface  Weapons  Center  (MSWC)  for  navagational  and  geodetic  satellites. 

The  program  accurately  determines  the  trajectory  of  a satellite  from  Doppler  beacon 
data.  CELEST  has  been  obtained  by  AFGL  and  modified  for  use  in  atmospheric  den- 
sity determination  Modifications  include  adaptations  to  the  operating  system 

(SCOPE  3.  4)  used  at  AFGL,  addition  of  the  Jacchia  1964^  and  other  density 
models,  and  the  determination  of  perigee  (position  where  geocentric  distance  mini- 
mizes) and  related  density  calculations. 

7.1  Results 

CELEST  was  used  to  analyze  data  from  satellites  DB-7,  DB-8,  and  DB-9  to  determine 

atmospheric  density  by  least  squares  adjustment  of  a scale  factor  (or  drag  coefficient) 

multiplying  the  model  density.  For  comparison  a parallel  study  was  performed  with 
(31 

program  CADNIP'  ' for  the  same  satellites  and  time  periods.  The  same  density 
model  (Jacchia  1964)  was  used  in  both  studies,  but  CADNIP  employed  a smaller 
geopotential  model  and  ADC  skin-track  observations. 

7.1.1  Perigee  Latitudes  and  Heights 

Consistent  disagreement  was  found  in  perigee  latitudes  and  heights  (Figures 
1-4)  for  satellites  DB-7  and  DB-8.  This  was  ultimately  traced  to  the  CADNIP  defi- 
nition of  perigee  as  the  location  of  zero  mean  mean  anomaly.  The  definition 
employed  in  CELEST,  given  previously,  is  equivalent  to  zero  osculating  -mean 
anomaly.  The  difference  between  osculating  and  mean  mean  anomaly  is 

= C?  X /8  ea2  (1  - e2) 
c x 

where  X = 2 (3  cos2i  -1)(1+X  ) (1-e2)  sin  E/(l-e  cos  E) 

+ 3 sfn21  [(l-)f  ) sin  (V  * 2W)  +{x ' + 1/3)  s1n(3V  + 2M) 

»1th  x'  - -hrcosT  t-i-'lasTE*1) 

E = mean  eccentric  anomaly 
V - mean  true  anomaly 
W = mean  argument  of  perigee 
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Figure  1.  Perigee  Latitudes  for  Satellite  DB-7.  Vertical  lines  show  thrust  times 


Figure  2.  Perigee  Heights  for  Satellite  DB-7 


Figure  3.  Perigee  Latitudes  for  Satellite  DB-8 


Figure  Perigee  Heights  for  Satellite  DB-8 


i = mean  inclination 
a = mean  semi-major  axis  (Earth  radii) 
e = mean  eccentricity 

* second  harmonic  coefficient 

O ' 

For  zero  mean  mean  anomaly,  values  of  the  other  elements  (W  = 135,  i = 96  , 
a = 6590  wn,  e = .009)  typical  for  these  DB  satellites , and  = -1082.645  x 107® 
there  results  4M  = 3.15°,  which  is  clearly  consistent  with  the  results  shown. 

CADNIP  was  therefore  modified  to  define  perigee  for  vanishing  osculating  mean 
anomaly.  The  improved  agreement  with  CELEST  (Figures  5-6)  for  satellite  DB-9 
verified  the  above  conclusion. 

7.1.2  Atmospheric  Density 

As  shown  in  Figure  7,  densities  obtained  from  the  two  programs  were  gen- 
erally in  good  agreement.  The  largest  differences  occur  near  times  of  thrusts,  the 
parameters  of  which  can  be  optimized  with  least  squares  adjustment  cy  CELEST 
but  not  by  CADNIP.  Correlation  of  density  model  ratio  increases  with  large  geo- 
magnetic activity  increases  (Figure  8)  are  seen  for  both  sets  of  results,  indi- 
cating a deficiency  in  the  density  model. 
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Figure  5-  Perigee  Latitudes  for  Satellite  DB-9 


Figure  7*  Model  Ratios  for  Satellite  DB-9 
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Section  8.  Geopotential  Model  Studies 
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8.0  Geopotential  Model  Studies 

Initiator:  K.  Champion 

Project  No.:  427M  Problem  No.:  4867 

A crucial  element  in  satellite  orbit  determination  is  the  selection  of  force 
models.  The  most  significant  forces  effecting  an  artificial  close-earth 
satellite  are  the  Earth's  gravitational  attraction  and  atmospheric  drag.  The 
determination  of  the  gravitational  field  (geopotential)  in  terms  of  spherical 
harmonics  is  straightforward;  however  a large  number  of  terms  are  required  for 
accurate  representation.  On  the  other  hand,  lack  of  knowledge  of  the  dynamics 
of  time  dependent  atmospheric  density  variations,  particularly  in  the  presence 
of  geomagnetic  disturbances,  has  been  a major  stumbling  block  in  accurate  pre- 
diction of  satellite  position. 

It  is  therefore  apparent  that,  in  many  cases,  particularly  low-altitude  satellites, 
the  use  of  large  geopotential  expansions  may  not  be  justified  in  view  of  com- 
putation time  and  memory  requirements,  save  for  the  inclusion  of  "resonance" 
terms.  It  is  therefore  desirable  to  know  the  acci  icy  of  various  published  geo- 
potential expansions  (models),  the  effects  of  tuncating  a model,  and  the  effects 
of. various  resonance  terms. 

The  work  reported  here  is  a continuation  of  the  geopotential  model  evaluation 
underway  at  AFGL  as  described  by  Garrett,  et  al.^  Satellites  DB-7  (6382), 

DB-8  (6727),  ini  DB~9 ( 6928)  have  been  chosen  for  this  study,  because  of  the 
availability  of  nigh  quality  doppler  beacon  data  at  AFGL.  Thus  a comparison 
of  results  for  different  types  of  data  is  possible.  Several  geopotential  models 
are  compared  as  to  their  orbit  determination  and  prediction  capabilities.  Modi- 
fications by  truncation  and  exclusion  of  resonance  terms  are  also  considered. 
Uncertainties  in  Earth  contents  (gaupctential  constant,  radius,  flattening  ratio) 
are  also  sources  of  errors  in  orbit  determination  and  prediction,  and  are  also 
explored  in  this  report. 

8. 1 Procedure 

Doppler  beacon  data  obtained  from  Naval  Surface  Weapons  Center  (NSWC),  Dahlgren 

(2  3) 

Laboratory  was  processed  by  program  CELEST'  ’ . Data  were  fit  over  approximately 

27  hour  time  spans  by  least-squares  adjustment  of  6 orbital  parameters  and  a drag 
factor.  A second  fit  was  obtained  in  each  case  for  a time-span  beginning  at 
approximately  the  midpoint  of  the  first  time  span  and  extending  to  16  hours 

after  the  end  of  the  first  span.  Prediction  errors  were  thus  obtained  from 
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positional  differences  between  the  orbit  obtained  from  the  first  time  span  and 
that  obtained  from  the  second. 

These  procedures  were  carried  out  for  the  geopotentials  listed  in  Table  1 and 

described  in  Reference  1,  and  for  time  spans,  selected  for  varying  levels 

of  geomagnetic  activity,  listed  in  Table  2.  Constants  were  input  as  prescribed 

for  each  model.  Tests  described  later  however  showed  that  variations  of  these 

parameters  over  the  ranges  of  published  values  produced  little  change  in  the 

(41 

results.  The  Jacchia  1964'  ' atmospheric  density  model  was  used. 


Table  1.  Geopotential  Models 

Model Maximum  Degree 


Gaposchkin 

8 

STEM  68R 

15 

SEIII 

23 

MGS  72 

28 

14th  Aerospace 


12 


Table  2.  Satellites  and  Fit  Spans  Studied 


Satellite  Avera9e  KP 


(perigee  Height) 

Fit  Span  (1973) 

Fit 

Prediction 

March  14 

1.0 

1.0 

DB7 

March  16 

2.5 

2.0 

(157-162  Km) 

March  18 

4.0 

6.0 

March  23 

5.5 

5.5 

July  26 

4.5 

4.0 

DB8 

August  3 

1.5 

2.0 

(162  Km) 

August  8 

2.5 

1.5 

August  10 

1.5 

1.0 

November  14 

2.0 

2.0 

DB9 

November  29 

1.5 

1.5 

( 1 64-1 68Km) 

December  6 

2.5 

2.0 

December  22 

4.5 

4.0 
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8.2  Results 


Prediction  errors  obtained  are  typified  by  Figures  1-3.  The  most  visible 
features  are  the  general  increasing  trend  of  the  total  and  in-track  errors, 
upon  which  are  superimposed  oscillations  with  a period  of  one  orbital  revolution. 
These  are  due  to  phase  and  amplitude  differences  in  the  oscillation  of  true 
anomoly  about  mean  anomoly  for  the  two  approximately  elliptical  orbits  being 
differenced  (See  Section  8.3).  They  are  removed  by  plotting  differences  in 
mean  anomoly  plus  argument  of  perigee  (Figure  3).  Errors  in  components  per- 
pendicular to  in-track  were  generally  less  than  200  m. 


8.2.1  Statistical  Analysis 

The  total  error  has  been  analyzed  as  described  in  Reference  1.  The 
following  quantities  have  been  computed  for  each  model  at  5,  10  and  15  hours 

beyond  the  fit  span 
— N 

ER  =1  I ERj 
N i«l  1 


SD  = 


M 


n — o 

I (ER.  - ER)2 

A =1  1 

N-l 


I er| 


lER1 1 


where  ER..  is  the  total  error  for  the  ith  case,  having  the  sign  of  the  in-track 
error. 

Tables  3-5  contain  the  results  for  12  cases  processed  by  program  CELEST.  As 
expected,  the  mean  total  error  ETT  is  generally  small.  The  more  significant  quantities^ 
the  standard  deviation  SD  and  the  mean  magnitude  of  error  |ER|  both  indicate  a 
marginal  preference  toward  the  higher  order  models  SE  III  and  WGS72.  This  is  in 
general  agreement  with  the  results  of  Reference  1,  except  for  the  good  performance 
of  the  SE  III  and  the  poor  performance  of  14th  Aerospace  12  x 12  model.  The 
errors  obtained  here  are  also  much  smaller  than  those  found  in  Reference  1. 
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ELAPSED  TIME  (HOURS) 


ELAPSED  TIME  (HOURS) 


This  could  be  due  to  use  of  more  accurate  data,  choice  of  longer  fit  spans,  or 
the  choice  of  time  periods  in  the  present  cases  for  which  atmospheric  drag 
error  was  relatively  small.  Nevertheless  in  inter-comparing  models,  we  reach 
the  same  conclusion,  that  for  drag  perturbed  orbits  the  small  improvement  in 
accuracy  gained  from  use  of  the  higher  order  models  is  outweighed  by  the  increase 
in  time  and  memory  requirements. 

Table  3.  Mean  Total  Prediction  Errors (km)  - Doppler  Beacon  Data 


STEM  68  R .241  .168  .257 


Gaposchkin  .600  1.395  1.332 


Table  5.  Mean  Absolute  Value  of  Prediction  Error  (km)  - Doppler  Beacon  Data 


: 


| 


K 


! 


Gp.  Model 

5 HR 

10  HR 

15  HR 

STEM68R 

1.214 

2.499 

4.193 

Gaposchkin 

1.181 

2.595 

4.436 

SEIII 

0.930 

1.822 

3.005 

WGS72 

1.004 

2.084 

3.699 

12  X 12 

1.266 

2.630 

4.599 

8.3  Comparison  of  Close  Orbit  Characteristics 

Studies  of  satellite  orbit  determination  programs  where  the  force  model,  or 
the  set  of  observations  and  the  time  span,  or  the  program  itself  is  varied 
result  in  close  but  not  identical  orbits  which  must  be  compared  to  evaluate  the 
effects  of  the  variations. 

Comparison  consists  of  determining  total,  in-track  (ITE),  radial,  and  cross-track 
errors  between  the  observed  (or  reference)  orbit  and  the  calculated  (or  predicted) 
orbit.  Since  the  trajectories  are  known  in  terms  of  geocentric  ECI  position 
and  velocity  as  a function  t"'me,  these  errors  may  be  directly  computed  as 
below. 


Let  the  observed  and  predicted  position  and  velocity  vector  components  be 
Xn  Y„  Zn  K Y L 

OOO  000 

Xn  Y 1 Xn  Y Zn 

P P P P P P 


Let  V, 


n~rr~ 

/ X 0 Yo  + 
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The  radius  vectors  are  given  by 


Then, 

Total  Error  = 


J 


In-Track  Error  including  sign  (ITE) 


J <xo  - V Xo  * <Yo  - V Yo  + <Zo  - V Zo 


Radial  Error  e R„  - R 
o p 


Cross-Track  Error  (CTE)  = R X + R , Y„  + R,  Z„ 

x o y o z o 

v 


Where  Rx  = YQ  Zp  - Y,  Zq, 
Ro 


Rv  = Zo  Xo  - Zo  V 
Ro 


(2) 

(3) 

(4) 


= Xn  Yn 

JL+ 


Positive  CTE  is  in  the  direction  of  the  orbit  spin  axis. 

Note  that  the  observed  and  calculated  velocity  components  may  be  used  interchange- 
ably as  available,  with  negligible  change  in  results. 

It  should  be  observed  that  the  three  error  components  above  do  not  form  an 
orthogonal  coordinate  system.  This  however  was  not  relevant  to  the  studies 
made  here.  In-track  error  in  the  velocity  direction  gives  a direct  timing 
error  which  is  useful. 
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Aerospace's  TRACE  orbit  determination  program  defines  In-track  error  as 


(X„  - »„)  T„  + (Y„  - y Ty  * (Zn  - Z„)  T, 

/t  Z + T >d  + T ' 1 
' x y z 


A 

1 


Where  n designates  "nominal"  orbit  and  o designates  observed  orbit,  and 
the  vector  Is  given  by 


r x C > x K 


The  difference  between  these  two  definitions  of  In-track  error  becomes 
significant  only  for  highly  eccentric  orbits  where  the  velocity  vector  Is  not 
normal  to  the  radius  vector. 

In  drag  dependent  orbit  studies,  the  total  error  turns  out  to  be  almost  com- 
pletely due  to  the  In-track  error,  while  the  radial  or  altitude  error  and  the 
crt>ss-track  error  remain  consistently  small.  This  In-track  error  will  gen- 
erally Include  four  components: 

1)  An  average  In-track  error  - offset 

2)  Oscillations  In  the  In-track  error 

3)  A linear  Increase  In  the  error 

4)  An  exponential  component  In  the  error 

Figures  1 and  2 show  typical  total  and  In-track  error  computer  plots  where  fit 
to  observed  data  has  been  compared  to  predicted  orbit  over  a 16  hour  period  for 
five  different  geopotentlal  models.  The  oscillations  are  clearly  evident, 
and  the  other  three  components  are  present  In  varying  degrees. 

Orbit  determination  programs  generally  also  provide  classical  orbital  elements 
at  one  or  more  epochs.  It  Is  possible  from  a comparison  of  these  elements  to 
quanltatlvely  determine  the  relative  behavior  of  close  orbits  over  a number  of 
revolutions.  The  four  components  of  In-track  error  are  primarily  due  to 


differences  in  the  semi-major  axis  (a),  eccentricity  (e),  argument  of  perigee 
(w),  mean  anomaly  (m),  and  the  drag  factor  (CD).  The  inclination  and  the 
longitude  of  the  ascending  node  are  factors  which  primarily  affect  cross-track 
error. 

In-track  error  at  any  time  is  given  by 

ITE  = R ( aw  + av)  (5) 

where  R = radius  vector 

v = true  anomaly. 

The  difference  in  true  anomaly  av  may  be  expressed  in  terms  of  a M and  a e using 
the  relationships 

M a E - e Sin  E (6) 

tan  v » /l+e  tan  E 

7 7 (7) 


where  E is  the  eccentric  anomaly. 


If  we  let 


we  obtain 


? 


= tan  E 

7 


and 


AV 


l 


e cos  E 


) aM  + j(1^)  St"  E 


(1-e  cos  E) 


-e 

+e 


±1 


Ae 


(8) 


“ ' + Tr> 


The  differences  aw,  aM,  and  Ae  coupled  with  the  radius  vector  R oscillating 
at  the  orbit  period  thus  produce  the  in-track  error  and  its  variations. 


8.3.1  Estimation  of  Typical  In-track  Errors  for  Low  Altitude  Satellites. 

Application  of  equations  5 and  8 for  an  orbit  with  a=1.3  and  ea.2  shows  that 
differences  Aw  or  AM  of  0.01,  or  Ae  of  0.0001  result  in  ITE  magnitudes  of 
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about  1 Km.  This  error  is  accompanied  by  oscillations  that  peak  near  180°, 

0®,  or  90*  true  anomaly.  Further,  unless  (Aw  +a  M)  is  nearly  zero,  an  average 
ITE  offset  Is  also  associated  with  these  oscillations. 

ITE  due  to  Argument  of  Perigee  - aw 


ITE  max  = a (l+e)A  w at  M=180 


ITE, 


min 


= a (l-e)  a w at  M=0 


(9) 


ITE  due  to  Mean  Anomaly  - aM 


ITE, 


max 


ITEmin 


a/lTe 

JpS 

a /l-e 

^i+e" 


A M at 


aM  at 


M =0 


M=180 


(10) 


ITE  due  to  Eccentricity  - Ae 


/ — 

ITE  Amplitude  ~ a(3+2e^  A e 


(ID 


o o 

Max  for  v s 90  ; min  for  v z -90 


Differences  Aa  and  aCq  result  respectively  in  linear  and  exponential  growth  In 
ITE.  The  latter  phenomenon  is  essentially  caused  by  a linear  growth  In  Aa  due 
to  differing  drag  or  atmospheric  density  models  for  the  two  orbits. 


ITE  due  to  Semi -Major  Axis  - Aa 

aITE  per  revolution  * - 31T  Aa.  (12) 

Thus  a 0.1  Km  difference  causes  a lKm  per  revolution  build  up  in  ITE. 
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ITE  due  to  Drag  Coefficient  - i Cp 

The  monotonic  decrease  in  'a'  and  the  orbital  period  is  due  to  the  drag. 

The  ratio  Acd/Cd  gives  the  portion  of  this  rate  of  change  of  'a*  that  will 

persist  and  show  up  as  a near  parabolic  increase  in  ITE.  If  f is  the 
fraction  decrease  per  day  in  the  semi-major  axis,  and  n is  the  orbital 
period  in  revolutions  per  day, 

ITE  value  attained  in  one  day 


= ACo 


3 IT  f n a 
' 2 


(13) 


Thus,  for  a =1.3,  n=l  1 . 5 revs/^  and  semi -major  axis  decay  rate  of 

0.1%  per  day,  a 1%  difference  in  the  drag  coefficient  will  show  up  as  an 
ITE  of  4.  5 Km  after  one  day. 


8.4  Resonant  Terms 

It  is  evident  from  the  results  obtained  here  that  one  may  indeed  often  reduce 
the  size  of  the  geopotential  model  without  adversely  affecting  prediction 
accuracy.  To  gauge  the  extent  to  which  a geopotential  can  be  reduced  and  tocsee 
which  terms  may  be  omitted  it  is  useful  to  have  a method  of  estimating  the  influence 
of  individual  geopotential  terms  on  prediction  accuracy. 

8.4.1  Geopotential 

As  stated  previously  the  geopotential  is  expressed  as  an  expansion  of  spher- 
ical harmonics 
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where 


GM  * Earth's  gravitational  constant 

r = radius  vector 

a^  = radius  of  earth 
e 

X * latitude 

Y = longitude  (East) 
p m 

= associated  Legendre  function 


8.4.2  Influence  of  Individual  Terms 

The  effects  of  individual  geopotential  terms  on  satellite  orbits  fall  into 
2 classes^ 

1)  Secular  (m=o  only) 

2)  Periodic 

A cos  2 TTf  t 

f = (j^2p)  W + U-2p+q)  n+m  ( n-  0) 

where 

„ = precession  in  argument  of  perigee 

n = precession  in  ascending  nod 
9 = sidereal  rotation  rate  of  Earth 

p = integer  (o  4 P 4 XI 

q -^integer  j£_ 


e = eccentricity 

The  most  important  of  the  secular  corrections  comes  from  theJ@s2,  m=o  term  since 
c20  is  1000  times  as  big  as  the  higher  order  coefficents.  These  include  the 

precession  of  perigee  and  ascending  node  appearing  in  the  periodic  effects. 


117 


The  amplitudes  A of  the  periodic  effects  are  functions  of  the  orbital  elements, 
which  are  slowly  varying;  hence  A may  be  considered  constant  over  sufficiently 
short  time  spans.  Since  it  and  {i  are  much  smaller  than  n or  e it  is  evident  that 
resonance  occurs  for 

n (revolutions/day)  = m. 

Thus  the  DB  satellites  (16.2  revs/day)  are  near  resonant  for  m=16  terms  in 
the  geopotential. 


8.4.3  Effects  on  Predicted  Orbit 

The  resulting  influence  of  a particular  geopotential  term  on  a predicted  orbit 
based  on  least-squares  adjustment  of  orbital  parameters  is  given  by 


= i**1  orbital  parameter 
C = coefficient  of  geopotential  term 


The  first  term  accounts  for  the  dependence  on  the  geopotential  of  the  orbital 
parameters  obtained  in  a least  squares  fit,  while  the  second  term  is  directly 
obtainable  from  the  perturbation  theory  described  in  8.4.2. 


8.4.4.  Dependence  of  Orbital  Parameters 

It  is  assumed  that  the  change  in  orbital  parameters  due  to  omission  of  a 
term  in  the  geopotential  will  be  such  as  to  provide  a lease  squares-adjustment 
to  the  perturbations  in  the  orbit  due  to  that  term. 

Since  the  dominant  prediction  errors  are  in-track  (Figures  1-3)  we  concen- 
trate on  this  component.  We  assume  a 3-parameter  fit  function  for  in-track 
position. 

3 

s/t)  = T p,  t1"1 
x=i  1 

with  P^representing  the  initial  orbital  longitude,  P?  proportional  to  the 
mean  motion  (function  of  the  semi-major  axis),  and  P3  proportional  to  the  rate  of 
change  of  mean  motion,  as  produced  by  atmospheric  drag. 

It  is  evident  that  secular  corrections  mentioned  in  8.4.2  will  be  fit  per- 
fectly in  so  much  as  they  are  limited  to  first  time  derivatives  of  the  elements. 
Their  effects  on  the  predicted  orbit  would  thus  cancel  out.  We  therefore  con- 
centrate on  the  periodic  corrections. 

Then  we  find  that  the  changes  dPi  due  to  a periodic  orbital  term 
A cos  (at  + qQ)  are 
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dP,  = - A 
J 3 

ar 

|30  (sin  q - 

sin  qj  + 180 
0 “if 

(cos  q + cos  qQ) 

where 

-360 

(aT)2 

(sin  q - sin  q 

T = length  of  fit  span. 

q = aT  + qQ. 


8.4.4  Results 

Figure  4 shows  the  rms  change  in  prediction  error  (averaged  over  phase  qQ). 

It  is  seen  to  depend  on  the  ratio  of  fit  span  T to  the  period  of  the  perturbation, 
and  the  ratio  of  the  prediction  time  measured  from  the  end  of  the  fit  span  to  the 
fit  span  T.  The  effect  on  prediction  error  is  largest  for  those  perturbations 
with  periods  nearly  equal  to  the  fit  span.  This  is  illustrated  in  Table  6. 
Satellite  DB7  has  a mean  motion  of  16.2  revs/day  and  is  thus  near-resonant  for 
m=16  terms  (principal  perturbation  has  5 day  period.).  The  m=17  terms  contribute 
a perturbation  with  a 30-hour  period. 
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Table  6.  Influence  of  m = 16  and  17  Terms 
on  DB-7  Prediction  Error 

Fit  Span:  27.5  hours  ending  at  0 hours  day  76  (1973) 
Predict  Time:  17  hours 
Geopotential:  SE  III 


Estimate  (Km) 

Actual  (Km) 

m 

(CELEST) 

16 

- .3 

- .5 

17 

-2.3 

-2.5 

8.5  Uncertainties  in  Earth  Constants 

Values  of  the  various  Earth  constants  specified  for  the  various  geopotential 
models  generally  differ.  Thus  it  would  seem  desireable  to  use  the  set  of  constants 
appropriate  to  the  particular  model.  Deviations  from  this  are  common  and  thus 
it  is  useful  to  ascertain  their  effects.  Table  7 lists  various  sets  of  values. 
Table  8 lists  16-hour  predicted  orbit  differences  obtained  from  program  CELEST. 

It  is  apparent  that  the  effects  of  varying  the  Earth  constants  over  the  un- 
certainties implied  in  Table  7 are  negligible  compared  to  the  total  errors  pre- 
sented earlier. 
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Table  7.  Earth  Constants 


aea  (km) 

l/fb 

GM 

(km3/sec2) 

Gaposchkin 

6378.145 

298.25 

398601.2 

STEM68R 

6378.145 

298.25 

398601.2 

SE  III 

6378.140 

298.256 

398601.3 

WG$72 

6378.135 

298.26 

398600.5 

14th  Aerospace  12  x 12 

6378.145 

298.25 

398601.2 

aEquatorial  Radius 
^Flattening  Ratio 

Table  8.  Eff: :t  of  Constant  Changes  on  SE 

III  Fit  and  Predict 

A a/  (m) 

5 

-5 

5 

Af _1 

- . f"  JC 

.004 

-.006 

Am  Uon?/sec2) 

0 

0 

- .8 

^x0b  ('«,.) 

1 

-1 

4 

A Yq  ( m) 

0 

1 

-3 

AZ0  ( DJ  •> 

-1 

-1 

-1 

AX,C  (.'  m.) 

2 

5 

1 
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Table  8.  (continued) 


AY1  ( m ) 0 -3  2 

AZ,  ( m ) -1  0 4 

a Measured  from  SE  III  constants  (row  3 of  Table  7) 
b From  fit  initial  conditions  8/11/73  20  hr  40  min.  for  satellite  DB-8 
c From  predicted  position  8/13/73,  16  hr.  for  satellite  DB-8 
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8.6. 


NORMALIZATION  OF  GEOPOTENTIAL  COEFF I C I ENTS 


Geopotential  coefficients  for  oroit  determination  programs  are  encountered 
with  3 types  of  normalization  in  common  practice: 


1)  Unnormalized 

2)  APL  or  partially  normalized 

3)  Kaula  or  fully  normalized 

Expressions  relating  the  coefficients  for  the  various  normalizations  are: 


Nf 


(n+mK 

(n-m): 

(n+m)! 


NJ  

(2n+l ) (2-6^)  (n-m)! 
1 

N (2n+l ) (2-6  ) nm’ 

MtO 


u 

nm. 


Where  the  superscripts  U,  A,  K imply  unnormalized,  APL-normalized  and  Kaula- 
normalized  coefficients  respectively,  and 

Vo 

* 0,  m f 0 

The  unnormalized  coefficients  only  are  required  by  CELEST  and  CADNIP,  while  any 

:\ 

of  the  3 may  be  used  in  TRACE  provided  the  flag  NFORM  is  properly  set.  The  APL 
normalization  was  used  by  Guier  and  hence  applies  to  the  Guier  deck. 

Most  labs  (SAO,  Spacetrack,  iowC;  -eem  to  prefer  the  Kaula  normalization.  Hence 
it  must  be  determined  what  normalization  has  been  used  for  each  deck  transmitted 
to  us  from  the  outside. 

A ■■•ord  of  warnino:  in  some  cases,  the  normalization  of  zonals  (m*0)  is  not  the 

same  as  non-zonals,  rather  the  J coefficients 

« r. 

J = -c u 
n no 

are  given  while  APL  or  Kaula  normalization  is  used  for  the  non-zonals.  Thus 
the  well-defined  (2,0)  coefficient  should  always  be  checked  in  addition  to  the 
normalization.  ,oc 
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Section  9.  Rocket  Trajectory  System 


t 


9.0  Rocket  Trajectory  System 


r 


Initiator:  E.  Robinson 

Project  No.:  0001  Problem  No.:  4517 

In  parallel  with  a continuing  active  use  of  the  AFGL  Rocket  Trajectory 
System,  requirements  for  improvements  in  nearly  all  phases  of  this  system 
were  identified.  A combined  report  and  user's  guide  is  in  preparation^ 
in  order  to  present  a current  document  of  procedures  and  pertinent  analytical 
backup.  Significant  additions  and  changes  since  the  last  report^  are 
outlined  here. 

9.1  Editor  - Preprocessor  (DRIVEA) 

Occasionally  the  raw  radar  data  contains  biases  or  extended  deviations  in 
range,  azimuth  or  elevation  which  are  discernible  by  Inspection  of  the 
preprocessed  plots.  A capability  was  implemented  in  DRIVEA  to  specify  one 
or  more  time  spans  over  which  it  is  desired  to  bridge  any  one  or  more  of 
these  measurements  by  a smooth  spanning  cubic  polynomial.  The  cubic 
satisfies  readings  at  the  two  times  specified  plus  two  readings  5 seconds 
outside  the  span  of  discarded  readings. 

9.2  Filter  (ORIVEB) 

The  main  developmental  effort  took  place  in  this  area  where  statistical 
filtering,  rocket  dynamics  modeling,  digital  integration,  and  output 
considerations  are  combined. 

9.2.1  Operating  Modes 

Figure  1 shows  the  types  of  operation  of  DRIVEB  that  are  available  by  input 
of  the  variable  KEYOPT.  Inputting  of  TB,  TC,  TE,  and  TF  relative  to  TLNCH 
and  TSTART  establishes  regions  for  the  consecutively  numbered  integrating, 
filtering,  or  reproducing  phases.  Geocentric  position  and  velocity  initial 
conditions  are  required,  except  that  launch  azimuth  and  elevation  can 
replace  velocity  information  at  launch.  KEYOPT  = 3 provides  a run  identi- 
cal to  the  original  DRIVE  3 run,  and  KEYOPT  = 4 matches  the  original 
DRIVE5.  KEYOPT  = 1 is  normally  used  now  and  is  discussed  below. 

9.2.2  Integration  - Filtering  Mode 

Reliable  radar  data  is  generally  not  available  till  some  time  into  launch. 

The  original  DRIVE  3^  and  DRIVE  5^  programs  filtered  data  from 

some  preselected  time  on,  and  attempted  to  obtain  the  trajectory  from  launch 

to  that  time  by  integrating  backwards. 
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This  procedure  inevitably  leads  to  an  unresol vable  mismatch  at  launch. 
The  solution  to  this  is  to  search  for  proper  launch  conditions  so  as 
to  achieve  the  closest  possible  match  to  subsequent  radar  position 
information.  Then  "pseudo"  radar  data  is  generated  from  launch 
allowing  a full  filtering  run  with  no  backward  integration.  The  search 
for  optimum  launch  determines  thrust  ratio  (fraction  of  nominal), 
launch  azimuth  and  elevation,  and  time  delay  before  onset  of  the  thrust 
phase.  The  square  sum  of  the  residuals  of  geocentric  coordinates  at 
three  time  points  near  the  start  of  good  radar  data  is  minimized 
iteratively. 


Briefly,  let  m be  the  number  of  variables  being  optimized  and  n the 
number  of  observations  (n>m).  Let  x0  be  the  vector  of  initial 
settings  of  the  m variables.  The  corresponding  residual  vector  for 

the  n observations  is  (rQj  r rQn).  m seParate  trials  are  carried 

out  with  small  increments  in  the  m variables  one  at  a time.  Let  these  increments 

be  given  by  Aq^  4qm,  where  q^  is  in  the  units  of  xi . The  augumented 

(m+l)xn  residual  matrix  R is 


11 


R = 


l 


Lrol 


In 


mn 


Solve  for  q from 


RR 


Lqnj 


= 0 


qm  1 


where  R is  the  transpose  of  R and  q0  + qj  + 
since  the  optimization  is  to  be  an  adjustment  about  x 

o 

Then,  if  the  are  changed  by  q^x..,  the  original  residual  vector  re 

will  be  minimized  in  a least  squares  sense.  Since  the  rocket  trajectory  is  a 
non-linear  process  the  optimization  is  generally  iterative.  The  m trials  pro- 
vide the  necessary  partial s without  analytic  methods,  and  are  considered  ap- 
plicable to  successive  iterations.  Taking  each  of  the  rectangular 
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coordinates  as  an  independent  observation,  the  result  is  the  minimization 
of  spatial  separation  between  the  observations  and  the  computed  trajectory. 

Highly  satisfactory  results  have  been  obtained  with  consistent  values  of 
thrust  ratio  and  delay  when  fitting  to  the  radar  data  at  different  times. 
Additional  process  "unknowns"  such  as  wind  can  be  added  to  the  search  as 
experience  identifies  and  warrants  these  variables. 

9.2.3  Driving  Noise  Uncertainty  £ 

The  driving  noise  uncertainty  in  the  velocity  vector  was  originally  taken 

2 

to  be  proportional  to  (ACCxDT)  where  ACC  is  the  instantaneous  acceleration 
and  DT  is  the  integration  step  size.  Corresponding  position  vector  un- 
certainty is 

ftp  v | (DT2/4  ) | 

It  was  found  that  the  dynamic  model  was  more  unreliable  at  re-entry  since  the 
integration  - filtering  mode  alleviated  the  launch  data  uncertainty. y is 
now  taken  as  proportional  to  W«t|ACC|»DT2  where  OVR  is  the  vehicle  range, 
ftp  is  calculated  as  before. 

9.2.4  Rocket  Dynamics  Modeling 

Drag  coefficient  determination  as  a function  of  Mach  number  has  been  refined. 
The  universal  curve 


CD=  l-0+'*mach 


V 


4.0 
^mach  + 


(Vmach^ 


^mach*1* 


is  multiplied  by  a rocket  and  stage-dependent  C0  multiplier.  Vmach  is 
given  by  vroC|<et/Vsound  ‘ ,J-S-  Standard  Atmosphere,  1962  provides 

a piece-wise  linear  estimate  of  the  speed  of  sound  as  a function  of  altitude 
and  this  has  been  incorporated  satisfactorily. 

Thrust,  mass  and  drag  characteristics  of  additional  rockets  were  added  to 
the  file.  The  repertory  with  rocket  key  numbers  is  given  in  Table  1.  Test 
integration  runs  and  adjustments  have  been  made  in  several  cases  based  on 
computer  runs  provided  with  the  nominal  specifications. 
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Rocket 


Key 


AEROBEE  150  03 

AEROBEE  170  04 
NIRO  07 
NIKE  TOMAHAWK  08 
UTE  TOMAHAWK  09 
PAIUTE  TOMAHAWK  10 
BLACK  BRANT  IV-B  16 
BLACK  BRANT  V -A  17 
BLACK  BRANT  V -B,C  18 
JAVELIN  19 
NIKE  JAVELIN  20 
NIKE  HYDAC  11 
SERGEANT-EXCEDE  25 
ASTROBEE  D 30 
AEROBEE  350  35 
CASTOR  40 
CASTOR- LANCE  41 


ARBITRARY  UNPOWERED 


NONE 


9.2.5  Integration  Step  Size 


r w 


Some  raw  radar  data  tapes  were  received  with  samples  at  1 second  rather  than 
the  usual  0.1  second  intervals.  It  was  found  that  if  the  integration  step 
size  was  allowed  to  increase  to  this  data  rate,  there  was  a significant 
build  up  of  error  in  the  4th  order  Runge-Kutta  integration  procedure. 

Step  size  less  than  0.1  second  resulted  in  negligible  change.  Subroutine 
RK  was  consequently  modified  to  loop  a selectable  number  of  times  with  a 
specified  step  size.  The  maximum  step  size  is  set  at  0.1  second. 

9.2.6  Output  Report  Data  File 

The  format  of  the  output  tape  was  substantially  revised  in  order  to 
a)  include  all  possible  trajectory  variables  that  were  desired  b) 
allow  REGEN  to  run  using  WRITER  with  no  further  computations.  The  original 
TAPE  8 and  TAPE  9 data  has  been  concentrated  onto  TAPE  4,  with  a header 
that  minimizes  card  input  with  REGEN  or  DRIVE  C (previously  DRIVE  4). 

Table  2 gives  the  format  . Subroutine  GENER  8 now  calculates  the  trajec- 
tory and  residuals  in  geocentric,  launcher,  and  radar  coordinate  systems, 
instead  of  only  the  ones  selected  for  printing.  The  resulting  TAPE  4 
may  be  used  by  WRITER  either  immediately  in  DRIVE  B or  DRIVE  C,  or 
subsequently  in  REGEN.  (See  Table  2.) 

9.3  Multi-Radar  Trajectory  (DRIVE  C) 

Program  DRIVE  C was  developed  as  a revision  of  program  DRIVE  4V  ' . The 
card  input  -'"uirements  are  greatly  reduced,  using  instead  the  header 
information  c.i  TAPE  4 generated  in  program  DRIVE  B (TABLE  2). 

The  primary  function  of  the  multi -radar  program  is  to  combine  the  results- 
of  two  separate  radar  observations  of  a given  rocket  flight  to  give  a 
composite  trajectory.  The  methods  used  to  accomplish  this  data  combination 
are  linear;  therefore,  an  arD itru.y  number  of  different  sets  of  radar  data 
pertaining  to  the  same  flight  may  be  combined  by  repeated  execution  of 
the  program.  The  grah  -,a  1,  printed  report  and  file  outputs  of  this 
program  are  identical  in  all  respects  to  the  outputs  of  Program  DRIVE  B. 

Specific  functions  of  Program  DRIVE  C include  the  following: 

• Process  two  input  files  to  produce  a single  composite 
trajectory  written  on  an  output  file 
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Table  2. 


Output  Report  data  file  written  on  TAPE  4.  The  first  record  contains 
flight  information  for  use  by  Program  Regen,  which  reproduces  printed 
and  plotted  output.  The  remaining  records  contain  data  to  be  printed 
and  plotted  by  either  DRIVE  C subroutines  or  by  REGEN.  If  I0PT=1,  all 
data  points  are  written  on  TAPE  4,  if  I0PT=0,  every  I SKI PTH  point  is 
written. 


FIRST  RECORD 


WORD(S) 

SYMBOL 

DESCRIPTION 

1-8 

TITLE  (1-8) 

80  character  alphanumeric 
identification 

9-10 

DAT (1-2) 

Run  date,  blank 

11-16 

TITLE  2(1-4) 

Further  identification 

TITLE  2(5-6) 

Payload  (KG),  Rarea  (KM**; 

17-20 

DAT(3) 

CD 

DAT(4) 

Launch  date 

21-23 

DAT ( 5-6) 

Rocket  ID 

24 

HRAD 

Height  of  radar 

25 

RLODEG 

26 

RLOMIN 

27 

RLOSEC 

Radar  position 

28 

RLADEG 

29 

RLAMIN 

30 

RLASEC 

31 

AALT 

Launcher  altitude 

32 

ALODEG 

33 

ALOMIN 

34 

ALOSEC 

35 

ALADEG 

Launcher  position 

36 

ALAMIN 

37 

ALASEC 

38 

UHRS 

39 

UMIN 

Universal  time 
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40 

USEC 

41 

DT. 

42 

TLI 

Initial  print  time  from  launch 

43 

TLF 

Final  print  time  from  launch 

44 

TSKIP 

Tape  written  every  TSKIP  secs. 

Table  2. 

(cont'd) 

SUBSEQUENT  RECORDS 

WORD(S) 

SYMBOL 

DESCRIPTION 

1 

GMT 

Universal  Time  (Hours) 

2 

TIME 

Time  After  Launch  (Seconds) 

3 

SN 

Signal  to  Noise  Ratio 

4 

RAG 

Right  Ascension  of  Greenwich(RAD) 

5-7 

P2( 1-3) 

Geocentric  Position  Components(Fi ltered) 

8-10 

P2 ( 4-6 ) 

Geocentric  Velocity  Components (Filtered) 

11-13 

ACC (1-3) 

Geocentri c Accel erati on ( Fi 1 tered) 

14-16 

ACC ( 4-6 ) 

Magnitudes  of  Geocentric  Vectors 

17-19 

0VR(l-3) 

Radar  Range,  AZ,  EL  (Filtered) 

20-22 

0VR(4-6) 

Radar  Range,  AZ,  EL  Rates( Filtered) 

23 

GEOALT 

Altitude  Above  Sub  Earth  Point 

24 

GEOLON 

Longitude  of  Sub  Earth  Point 

25 

GEOLAT 

Latitude  of  Sub  Earth  Point 

26-28 

GPV ( 4-6 ) 

Geocentric  Spherical  Coordinates 

29-64 

XSI ( 1-36) 

Complete  Error  Covariance  Matrix 

65-67 

DVR( 1-3) 

Raw  Data  Vector  Range,  AZ,  EL 

68 

RANGE 

Ground  Range  Along  Spheriod 

69 

RESALT 

Altitude  Residual  (Raw-Filtered) 

70-72 

XR.YR.z.P. 

Radar  Centered  Coordinates 

73-75 

XL,YL ,ZL 

Launcher  Centered  Coordinates 

76 

RLMAG 

Filtered  Launcher  Vector  Residuals 

77-79 

SRSG,YRSG,ZRSG 

Geocentric  Position  Vector  Residuals 

80 

GRSG 

Geocentric  Vector  Magnitude  Residual 

81-83 

XRSL,YRSL,ZRSL 

Launcher  Position  Vector  Residuals 

84 

GRSL 

Launcher  Vector  Magnitude  Residual 

85-87 

XRSG.YRSG.ZRSG 

Radar  Position  Vector  Residuals 

88 

GRSR 

Radar  Vector  Magnitude  Residual 

89-91 

P2 ( 4-6 ) 

Geocentric  Velocity  Component 
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Table  2.  (cont'd) 


Launcher  Velocity  Conponents 


92-94 

DVG(4-6) 

95-97 

OVR ( 4-6 ) 

98 

RRMAG 

90 

ALT 

*> 
1 
■ * 


Radar  Velocity  Components 
Radar  Range  Magnitude  |=0VR(  1 )J 
Raw  Altitude 


— Determine  the  type  of  data  region  being  read  for  each  file 

— Use  the  region  types  to  select  proper  mode  for  combination 
equations 

— Evaluate  combination  equations,  depending  on  the  particular 
mode,  to  produce  the  composite  observation  vector  and 
covariance  matrix 

--  Complete  the  calculations  for  the  remainder  of  each  output 
file  record  using  composite  values 

— Generate  a printed  summary  of  the  combination  process,  i.e., 
an  event  log  for  the  resulting  output  file 

• Produce  an  output  tape  file  containing  the  composite  trajectory 

• Generate  printed  and  graphical  reports  in  the  same  format  as 
DRIVED. 

For  each  radar  tracking  coverage  of  a rocket  a conventional  DRIVEB  tra- 
jectory vector  P(k)  with  its  associated  error  covariance  matrix  t(k)  is 
generated  for  each  time  point  k.  Two  radar  tracks  of  the  same  rocket 
may  be  statistically  merged  using  DRIVEC  to  provide  the  best  estimate 

P(k)  -[^(k)  Pj(k)  +*21(k)  P2<k)] [♦j1(k)  ♦♦^(k)]"1. 

The  associated  joint  error  covariance  matrix  is  given  by 

* (k)  -[♦'I(k) 

This  program  now  operates  in  only  one  of  two  basic  modes.  When  filtered 
or  integrity  data  is  available  on  both  tapes  the  above  merging  procedure 
applies.  Filtering  and  integration  regions  are  not  distinguished  except 
through  the  error  covariance  matrices.  When  data  is  absent  on  one  tape 
the  parameters  from  the  other  tape  are  directly  copied.  Check-out  of  the 
system  was  underway,  promising  early  results. 


Operating  System 


The  rocket  trajectory  operating  system  was  modified  to  reduce  the  need 
for  complete  binary  programs  of  DRIVEA,  DRIVEB,  DRIVEC,  and  REGEN.  The 
SCOPE  - EDITLIB  facility  is  used  to  provide  one  binary  library  for  all 
subroutines.  The  main  program  with  the  LIBRARY  instruction  then  causes 
unsatisfied  externals  to  be  loaded.  Program  maintenance  thus  requires 
only  one  file  of  source  routines  in  UPDATE  form. 
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Section  10.0  Satellite  Wire-Boom  Dynamics 


Initiator:  M.  Smiddy 

Project  No.:  8617  Problem  No.:  4702 

For  the  purpose  of  control  information,  mode  frequencies  determination,  stability 

analysis,  and  damping  characteristics,  a comprehensive  mathematical  analysis  and 

digital  computer  simulation  on  coupled  satellite-boom  dynamics  has  been  performed. 

(1  2 3 4) 

Since  the  full  details  are  reported  elsewhere  v * * ’ , only  the  salient 

features  of  this  work  are  outlined  in  this  chapter. 

(5) 

The  dynamics  of  a satellite  boom  system  is  generated  from  its  Lagrangian'  . 

For  dynamics  in  the  spin  plane,  the  Lagrangian  L is  derived  as 

4 

L * L0  («  , X,  Y)  +E  L.  (0  r,,  r,,  fl.,  |L  ) + 
i =1 

4 i 

H L.  (jfc  r,,  r *• , 0,,  0.,  X,  Y ) , (1) 

i=l  1 iiii 


where  L0  is  the  Lagrangian  of  the  hub,  L.  is  that  due  to  pure  rotation,  is  that 

due  to  cross  effect  of  rotation-translation,  the  arguments  of  the  Lagrangians  are 
the  canonical  variables,  dots  denote  time  derivatives,  and  subscript  i labels 
boom  number.  The  time  development  of  the  dynamics  of  the  system  is  determined  by 
the  Lagrangian  equations  of  motion,  viz.. 


(a)  the  hub  equation:  d j£L 

dt 

- ^ = _ K Q 

J0  K0 

(2) 

(b)  the  boom  equations:  d 31 

dt  9$ 

- ji  • -yi  <i-i.  ...♦> 

(3) 

(c)  the  tension  equations: 

_d  a L 
dt 

- = - Ti  (i=1 » • • *4) 

ri 

(4) 

(d)  the  translation  equations: 

A 

dt  55 

. = 0 

(5) 
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and 


(6) 


d 

at 


ak 


dL 

dY 


= 0 


where  L is  given  in  Equ.  (1). 

The  algebraic  details  of  the  above  equations  are  rather  lengthy^ . However, 
these  coupled  nonlinear  differential  equations  can  be  solved  simultaneously  and 
numerically  using  Hamming's  predictor-corrector  method^'  The  solutions 
represent  the  time  development  of  the  canonical  variables  of  the  system,  depending 
on  various  given  initial  and  boundary  conditions. 

Figure  1 displays  the  development  in  time  of  the  canonical  variables  under 
typical  initial  and  boundary  conditions. 

Analytically,  the  Lagrangian  equations  Equs.  (2)  to  (6)  can  be  solved  for  harmonic 
oscillation  modes.  Matrix  formulation  of  the  eigenvalue  equation  is  used. 
Simultaneous  diagonal ization  of  the  kinetic  and  potential  energy  matrices  (with 
dimensions  7x7)  yield  the  eigenvectors  and  the  mode  frequencies.  Table  1 summarizes 
the  analytical  results  of  harmonic  (normal)  modes  in  the  spin  plane. 

Due  to  the  presence  of  a very  effective  wobble  damper,  oscillations  out  of  the  spin 
plane  are  negligible  However,  for  the  benefit  of  future  design  and  for  the 

worst  event  of  tamper  failure  in  flight  time,  out  of  spin  plane  dynamics  have 
also  been  analyzed.  The  methods  as  outlined  in  the  preceding  paragraphs  on  spin 
plane  dynamics  are  used.  Table  i summaries  the  results  of  harmonic  mode  out  of 
spin  plane.  Figure  2 shows  the  frequencies  of  various  in  plane  and  out  of  plane 
modes  as  functions  of  boom  length. 

For  each  mode  excited  at  an  initial  time,  the  transient  response  of  the  system 
can  be  obtained  by  using  Laplace  transforms.  The  forced  oscillation  is 

0]  = [f ( t ) 0(t)]  (7) 

where  D is  a nonlinear  time  differential  operator  obtained  from  the  Lagrangian,j£is 
the  Laplace  operator,  F(t)  is  the  forcing  function,  and  0(t)  is  the  step  function 
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ANGULAR  VELOCITY  OF  HUB 


Figure  1.  Simulation  of  Hub  Angular  Velocity  and  Boom  Deflections  Following  a Disturbance 
Extension  of  booms  (1,3)  and  simultaneous  retraction  of  booms’  (2,4)  by  the  same 


ble  1.  Normal  Modes  of  Satellite  with  Pairwise  Equal  Length  Wire  Booms 


Pur*  Translation  Mod*  (Trivial) 


Table  2.  Normal  Mode*  Out-of-spin-plane  of  Satellite  with  Four  Equal  Length  Wire  Booms  - 


Mr*  Tr*a*Utl*a  Mod*  (Trivial) 


defining  the  duration  of  the  exciting  force,  such  as  boom  extension/retraction. 
Insight  on  stability  and  damping  characteristics  can  be  obtained.  Figure  1 
shows  the  excitation  and  decay  of  the  boom  oscillations  due  to  200  sec.  of  boom 
deployment. 
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11.0  Ionospheric  Field  Mapping 


Initiator:  M.  Smiddy 

Project  No.:  8617  Problem  No.:  4702 

This  work  supports  AFGL  research  projects  on  satellite  sensing  and  mapping  of  the 
global  distribution  of  electric  fields  in  magnetospheric  regions.  Extensive 
analytical  and  computational  procedures  are  required,  and  this  work  is  still  in 
process.  Schematically,  the  scope  of  this  work  is  outlined  in  the  flow  chart 
in  Figure  1. 

Data  reduction  considerations  including  handling  of  various  instrumental  and 
environmental  factors  are  discussed  below.  Preliminary  graphical  presentations 
of  the  results  are  also  included. 

11.1  Data  Reduction 

For  data  reduction  and  computation,  three  types  of  input  data  are  required. 

They  are  1)  electric  field  experiment  data,  2)  satellite  sensor  attitude 

(0M)  data,  and  3)  ephemeri s/magnetic  field  (ORMAG)  data.  The  structures  of 

(1  2 31 

these  data  are  tabulated  in  the  references'  ’ ’ 

The  electric  field  experiment  data  are  in  the  fifth  file  of  every  tape  of  the 
satellite's  raw  data.  The  experiment  data  contain  six  channels  of  sensor  signals 
with  various  high  and  low  gains.  These  signals  are  given  in  telemetry  counts, 
which  can  be  converted  to  voltages  by  a simple  conversion  factor. 

Calibration  sequences  show  up  every  512  sec  to  enable  the  updating  of  temperature 
correction  coefficients,  vehicle  potentials,  and  plasma  sheath  impedance, 
and  the  checking  of  their  statistical  deviations  for  various  channels. 


11.2  Instrumental  and  Environmental  Factors 

Before  the  genuine  electric  field  vectors  along  the  trajectory  of  the  satellite 
can  be  computed,  various  instrumental  and  environmental  factors  must  be  recognized, 
treated,  and  processed.  Intermediate  variables  with  physical  significance  are 
thereby  developed  and  made  available  to  the  researchers. 
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Prefliqht  Test  Data  Post-launch  Data  Analysis  Theory  and  Simulation 

. Design  logical  programfor  . Statistical  scheme  to  handle  signal  . Lagrangian  formulation  of 

satellite  data  reduction  saturation.  satellite  boom  dynamics, 

and  analysis. 
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Figure  1.  Task  Outline  for  Ionospheric  Field  Mapping  Project 


11.2.1  Saturation  Signals 


A statistical  method  is  used  to  treat  saturation  signals,  which  often 
appear  after  a calibration  sequence.  This  method  provides  modified  results  on 
plasma  impedance, vehicle  potential,  and  the  first  zero  crossing  time  TQ  of  the 
sinusoidal  signals  following  the  end  of  a calibration  sequence. 

The  extrema  locations  T^  of  the  saturation  signal  period  can  be  determined  by 
using  a least  square  sine  and  cosine  fit  of  the  form: 

Y = a-j  + (a2  + Agt)  Sin  wt  + (a^  + A^t)  cos  ut 
The  least  square  solutions  are  obtained  by  solving  the  matrix  equation:  - 
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Where  s^  = sin  (wt. ) 
Cj  = cos  (wt^) 


The  extrema  T.  obtained  from  equation  (T)  are  used  to  fit  the  sine  equation: 

J 

V1(t)  = 2[c.  + di  (t  - TQ) J sin  ut  (t  - TQ)  + e.  ( 2 ) 


by  solving  the  matrix  equation: 


f'W  sj  Vj 


5(trTo)sj  5<y  V2  sj  f'yv  sj 

5sj  j(ti-To)sj  J 1 


Vi(tj)sj 


where  Sj  = sin  £w(tj-  TQ)j 


and  Tq  is  the  time  of  the  first  zero  crossing  given  by  equation  (2)  with  offset  ei 
removed,  and  is  given  by 


where  K='3  if  the  corresponds  to  a minimum,  and  K=1  if  T-|  corresponds  to  a 
maximum. 


The  plasma  impedance  ^R1/  is  calculated  as 
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and  j corresponds  to  the  jth  interval  of  a calibration  sequence.  Similarly, 
the  vehicle  potential  ^JO1  > is  given  by 

Ki<>  1.2=3 

»here  *1.1  ' ijVM-V1.2+  Vk,3-  Vi,4  - Ak,l+  M.2  - Ak,3+  <4 
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The  statistical  deviations  of  these  quantities  are  also  calculated. 


11.2.2  Effect  of  Earth's  Rotation 

The  electric  field  induced  in  the  dipole  sensors  of  a satellite  is  of  sine  wave 
form  because  of  the  spin  of  the  satellite.  In  the  regions  of  low  latitude 
(equatorial)  and  low  altitude,  it  is  well  known  that  electric  field  is 
small  due  to  lack  of  significant  plasma  motion,  so  that  the  predominant 
electric  field  is  due  to  the  Earth's  rotation.  The  two  laws  invoked  in  this 
case  are: 


1.  The  electric  field  £ measured  in  a moving  frame  is  related  to 
the  electric  field  in  a rest  frame  by^ 

£ = E*  + V x B (3) 

where  il  is  magnetic  field  in  the  rest  frame,  and  velocity  is  assumed 
small  (V«C)  otherwise  there  is  an  additional  relativistic  factor  in  the 
denominator. 

2.  The  velocity  V in  a rotating  frame  r is  related  to  the 
velocity  V in  a fixed  frame  byv  ' 

^ = V-Wxr  (4) 

where  u>  is  the  rotation  angular  velocity  and  r is  the  position  vector. 

Thus,  the  induced  electric  field  due  to  the  Earth's  rotation  is 

£ = ^ x B (5) 

where  V is  given  by  eq.(4). 

This  induced  field  £ is  always  subtracted  out  in  order  to  bring  out  the 
genuine  electric  field  from  the  processed  result. 

11.2.3  The  Case  of  the  Missing  Dipole 

Satellite  signals  show  that  one  pair  of  dipole  sensorson  the  satellite  S3-2 
is  not  functioning  properly.  This  poses  a problem  on  the  determination  of 
three  dimensional  electric  field  vector  because  one  component  is  always 
missing.  However,  since  the  satellite  is  spinning,  a method  can  be  worked 
out  to  overcome  this  problem.  When  the  electric  field  is  maximum  in  the 
xy- plane  defined  by  two  dipole  pairs  and  r2  on  a spinning  satellite, 
it  is  reasonable  to  assume  that  the  electric  field  in  the  z -direction  is  zero. 
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However,  rj  and  need  not  be  orthogonal  at  flight  time,  causing  the 
vector  algebra  manipulation  to  be  nontrivial.  The  result  is 

E = Eji  * { E;  - (l!.r  2)  A + 

1 * r 2) 

where  ^denotes  the  unit  vector. 

This  gives  the  electric  field  vector,  inspite  of  the  missing  dipole. 

11.2.4  Effective  Frequency  and  Boom  Lengths 

The  information  on  the  frequency  u> of  the  spin  of  a satellite  is  given 
by  Boeing's  OM  data  tape.  The  electric  field  measured  by  the  satellite 
dipole  sensors  is  predominantly  (VxB.)  where  the  velocity  V has 
to  be  corrected  for  the  effects  of  Earth's  rotation.  The  value  of 
frequency  uj  from  OM  data  often  turns  out  to  be  slightly  different  from 
the  electric  field  variation.  Hence,  a curve  fitting  is  carried  out 
in  case  one  needs  very  accurate  indication  of  u> . 

Furthermore,  the  amplitude  of  the  induced  electric  field  in  a dipole 
depends  on  the  attitude  and  boom  length  of  the  dipole.  However,  the 
information  on  these  items  as  given  by  the  satellite  signal  data  tape 
may  not  be  perfectly  accurate,  due  to  the  possibility  that  the  booms 
of  a dipole  may  not  be  colinear  and  that  the  spin  axis  tilts  and 
precesses.  These  uncertainties  of  satellite  information  have  to  be 
treated  by  curve  fitting. 

For  orbit  nV,  an  analytical  computation  has  been  performed  assuming 
that  the  non-col  inear  booms  of  a dipole  are  simply  perpendicular  to 
the  instantaneous  spin  axis,  with  offset  angles  as  supplied  by  Aero- 
space Inc.  The  effective  boom  lengths  and  attitude  agree  quite  well 
with  the  least  square  curve  fitting  of  the  function  R.,(VxBj  against 
satellite  data. 

11.2.5  DC  Offsets 

Analysis  of  the  data  of  satellite  S3-2  shows  considerable  DC  offsets 
in  the  electric  field  data.  This  phenomenon  has  never.been  noticed 
before  because  no  experiment  using  such  a high  sensitivity  instrument  had 
been  performed  previously.  Further  analysis  shows  that  there  exist 
correlations  between  DC  offsets,  electron  density  and  oxygen  ion 
population  in  the  satellite  environment.  Since  the  ambient  electric 

1 
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field  is  mainly  associated  with  drifting  plasma^  in  the  polar  and 
auroral  regions,  the  DC  offset  is  fitted  in  the  equatorial  region  by 
means  of  a high  order  polynomical  (e.g.  13th) order,  which  is  then 
subtracted  from  the  computed  E-field  results  which  has  Induced  field 
( V- Vx  r)x  B already  removed. 

11.3  Graphical  Presentations 

11.3.1  E-Field  Strip  Chart 

Figure  2 shows  the  computed  ambient  electric  field  during  an  orbit. 

A number  of  estimates  are  made  during  each  satellite  spin,  resulting 
in  the  bar  pattern.  The  total  and  components  in  ECI,  geomagnetic 
(GM),  and  solar  magnetic  (SM)  coordinates  are  plotted  together  with 
environmental  information. 

11.3.2  E-Field  and  Drift  Current  Polar  Maps 

A typical  mapping  plot  of  electric  field  vectors  in  the  polar  region 
. 2 

is  shown  in  Fig.  3,  and  the  corresponding  (ExB)/B  drift  current  pattern 
is  shown  in  Fig.  4. 

11.3.3  Combination  Strip  Chart 

It  is  quite  often  useful,  when  analyzing  data,  to  have  a graphic 
presentation  of  other  physically  related  parameters.  In  the 
E-field  experiment,  parameters  such  as  the  ambient  plasma  density, 
the  heavy-ion  density,  and  the  high-energy  (1-5  kev)  particle 
flux  are  of  particular  interest,  since  each  can  affect  the  measurement 
in  its  own  particular  way.  Another  parameter,  the  orientation  of  the 
measuring  antenna  w.r.t.  the  ambient  magnetic  field,  is  useful  in 
determining  whether  or  not  a particular  anomaly  is  due  to  particle 
streaming  along  the  field  lines.  A program,  BILD3,  has  been  developed 
to  merge  and  present  all  parameters  and  data  of  interest. 

Functional  npcrrintinir 

The  plasma  and  ion  densities  and  particle  flux  can  be  obtained  from 
on-board  experiments.  The  antenna-field  orientation,  however,  must 
be  derived  from  model-field  calculations  and  satellite  attitude  data. 
BILD3  accesses  the  pertinent  data,  ephemeris,  and  altitude  files, 
unpacks  the  data,  performs  the  necessary  conversions  and  tranformationsj 
plots  on  one  graph  the  above-mentioned  parameters  together  with 


Figure  2.  Computed  Ambient  Electric  Field  During  an  Orbit 


E-field  data,  and  labels  the  X-axis  with  geophysical  coordinates 
such  as  universal  time,  altitude,  local  time,  and  invariant  latitude. 
Although  short  stretches  of  time  can  be  plotted,  as  in  Figure  5, 
the  program  is  more  often  used  for  entire  orbits,  producing  graphs 
six  to  eight  feet  long. 
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12.0  Plasma  Bulk  Motion 


Initiator:  P.  Wlldman 

Project  No.:  8617  Problem  No.:  4701 

Disturbances  and  Irregularities  in  the  ionosphere  and  magnetosphere, 
particularly  In  the  polar  and  auroral  regions  affect  electromagnetic 
wave  communications,  detection,  tracking,  guidance,  and  early- 
warning  defense  systems.  Ambient  plasma  bulk  motion  Is  an  Important 
factor  In  forecasting  Ionospheric  and  magnetospherlc  plasma  meteorology 
and  In  understanding  the  processes  leading  to  polar  region  Irregulari- 
ties. This  work  supports  research  projects  at  AFGL  on  In  situ  satellite 
sensing  of  plasma  bulk  motion.  Schematically,  the  tasks  performed  for 
this  problem  are  outlined  In  Figure  1. 

12.1  Data  Processing 

For  data  processing,  three  data  tapes  are  required,  viz.,  the  plasma 
experiment  data  tape,  the  satellite  attitude  (OM)  data  tape,  and  the 
satellite  ephemerls  (ORMAG)  data  tape.  The  structures  of  these  data 
are  explained  In  the  references.^1,2) 

The  satellite  signals  consist  of  several  patterns  which  have  to  be 
recognized.  They  are  the  TM  applied  voltage  sweeps,  electron  sensor 
real  sweeps,  electron  calibrates,  Ion  calibrates,  and  range  switchings. 

The  TM  sweeos  are  doubly  filtered  to  remove  any  nolsydata  that  may 
be  present.  The  processed  data  are  used  In  the  computations  of  electron 
and  Ion  properties  In  the  Ionosphere  and  magnetosphere.  The  flow  chert 
of  Figure  2 and  the  supplemental  notes  present  the  data  processing 
and  computation  logic. 

12 . 2 Electron  Properties  >•*  sm 

A TM  applied  voltage  sweep  appears  once  every  128  sec.  In  the  experiment. 
This  voltage  is  also  applied  to  every  real  sweep  of  the  electron  sensor 
In  the  following  128  sec,  period.  The  voltage-current  response  charac- 
teristics of  the  sensor  output  gives  Information  on  the  electron  proper- 
ties of  the  plasma  environment.  Physically,  the  electron  sensor  current 
I Is  proportional  to  the  velocity,  thermal  distribution,  and  density 
of  the  electrons,  l.e., 

I*NeA  Jit!  exp  (-$) 
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Summary  of  TM  sweeps,  real  sweeps,  electron 
and  ion  calibrates,  electron  densities 
and  temperatures,  plasma  flow  results. 
Create  virtual  tapes. 

Graphical  presentations. 
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.Generate  trial  values  for  FL'EPOW. 


ie(=ig 


Nef.Me 


SENSOR 

.Average  electron  current  before  and 
after  a real  sweep. 

.Correct  electron  density  with  Flepow. 


NerNeffe 


CRAFT 

CURRENT 

iSMil 

.Transform  ECI  coords, 
to  spacecraft  coords. 

i-on  currents  of* 8 sensors. 
.Ave.  iuu  ^urr=nts. (Issec.  ave.) 
.Select  forward  sensor  group 
.Compute  In/I,  state.  , , 

.Compute  velocity  angles  VOX, VOX, VOX 

1 

B 

f ell  • 

[.Compute  angles 

1 . 

FLOW 

.Compute  plasma  flow  angles  in 

spacecraft  coords. 

SPEED 

.Compute  plasma  flow  velocity. 

Figure  2.  Data  Processing  and  Computation  Flow  Chart  for  Plasma 
Flow  Parameters 
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Supplement  to  the  "Plasma  Flow"  Flow  Chart 


The  OM  attitude  data  tape  gives  interpolated  data  at  any  required  time 
within  a time  span  (approximately  equal  to  the  orbit  period).  The 
following  relevant  variables  are  available  in  Earth  Centered  Inertial 
System  (ECI). 

P— Position  vector  (Nautical  Miles) 

V. Velocity  vector  (Nautical  Miles/Sec.) 

S Direction  cosines  of  sensor  L.O.S.  in 

ECI  coord,  etc.  Ref.:  Boeing  OM  Manual  (1)J  . 

The  ORMAG  ephemeris  data  tape  gives  interpolated  data  in  the  same 
way  as  OM.  The  following  relevant  variables  are  available. 

P_ Position  vector  (km) 

V. Velocity  vector  (km/sec) 

B Magnetic  field  (Gamma) 

L Geomagnetic  Coodinate  L. 

Sun Solar  zenith  angle,  solar  longitude,  solar 

right  ascension  & solar  declination  . 

M.L.T.,  L.T.,  LAT,  LONG,  INV.  LAT.,  etc. 


[Ref.:  Mclnerney's  AFCRL  report  ^ ' 


• Instrument  temperatures  T]  to  T^  are  defined  as:  Tj  for  sensor 
package  219-1,  T^  electron  amplifier,  T-j  sensor  package  219-2, T4 
electronic  package 

(Tj,  Tz,  T3)  » (VT  - 1. 680 )/0. 0448 
T4  (Vt  - 2.190)/0.0404 

Ref.:  Peter  Wildman's  notes,  March  1975 

• On  board  ratios  R.  . are  defined  as  (I  .-I .)/( I .+1 .) 

' J J ' J 

• The  applied  voltage  V during  a TM  sweep  ranges  from  8.00 

ap 

volt  to  -8.50  v.;  it  is  converted  from  its  analog  linearly, 
which  ranges  from  5.10  volt  to  -O.lOv. 

• Electron  current  I is  obtained  as  log  I = log  I + kV 
where  V is  the  data  of  electron  sensor.  Putting  in  constants, 
we  have 

Ie  = exp  | -23.8020759  + 0.0041976  T?  + [2.4691428 
+ -.0017304  T2J  .vj 
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Supplement  to  the  “Plasma  Flow11  Chart  (Conti nued) 

Electron  density  Ng  is  given  by 

Ne  = 4.3535  x 1011  Ig  / f(v) 

where  f(v)  is  chosen  as  1.00,  Tfi  as  2500. 

• Flepow  gives  corrected  electron  density  Ne^,  vehicle  potential, 
and  electron  temperature  Tg  after  every  real  sweep.  If  the 
electron  sensor  is  partially  saturated,  Flepow  rejects  the 
sweep. 

• Ig  is  average  electron  current  before  and  after  a real  sweep. 


(15.25  15.78)sec. , (18.75  19.25)  sec.,  if  Hi  rate; 
(31.25  31.75)sec.  (34.75  35.26)  sec.,  if  Lo  rate. 
^Ref.:  Wildman  notes.  Sept.  1975^ 


•Corrected  electron  density  Ngc  = 
•Ion  currents  Ij  to  Ig  are  given  by 


Ie 

Te 


i.  = jk.v  + c]  x io"M 

where  is  temperature  dependent  function,  tabulated  in 
Wildman' s notes,  July  1975;  V is  ion  sensor  data,  C is  chosen 
as  -0.1,  and  M is  a range  dependent  number. 

•Time  average  of  *on  current  is  optional,  from  1/16  sec.  to 
1 sec.;  currently  h sec.  average  is  used. 

• oJier  items  in  the  last  boxes  of  the  flow  chart  are  less  basic. 


details  of  these  items  are  referred  to  Wildman 's 
write  up  on  April,  1976 
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where  the  square  root  term  is  due  to  the  electron  velocity,  the  ex- 
ponential term  is  the  Boltzmann  distribution  function,  N is  electron 
density,  e is  charge,  and  A is  sensor  aperture  area.  Therefore,  the 
slope  of  loge(I)  plotted  versus  voltage  V(t)  should  yield  the  electron 
temperature  T,  and  once  T is  found,  N can  be  computed.  However, 
in  reality  it  turns  out  that  the  current  response  of  the  instrument 
Is  highly  sensitive  only  in  a very  short  time  span  near  the  middle 
of  the  sweep  period.  Thus,  a slope  measurement  becomes  non-trivial. 

In  order  to  perform  highly  accurate  computation,  a subroutine 
FLEPOVI  is  used,  employing  the  steepest  gradient  minimization  technique 
of  Fletcher  and  Powell It  yields  results  on  electron  density, 
electron  temperature,  vehicle  potential,  and  plasma  sheath  radius, 
for  every  sweep.  Calculated  electron  properties  for  a typical  polar 
orbit  are  shown  in  Figure  3. 

12.3  Ion  Properties  and  Flow 

Eight  sensors,  divided  into  two  groups,  provide  measurement  data  of 
ion  currents.  It  is  the  forward  looking  group  which  gives  more 
accurate  information.  Since  the  satellite  is  always  spinning  the 
two  sensor  groups  alternate  in  facing  forward.  Thus,  sensor 
attitudes  and  velocity  vectors  are  computed  at  all  times  to  determine 
which  sensor  group  to  be  processed  at  every  moment.  For  diagnostic 
purpose  ground  computation  of  ion  sensor  current  ratios  are  compared 
with  on  board  ratios.  To  measure  bulk  ion  plasma  flow,  there  are 
approximately  two  regimes,  viz,  the  relatively  fast  ion  regime  and  the 
relatively  slow  ion  regime.  The  former  refers  to  the  light  ion  mass 
regions  such  as  those  occupied  predominantly  by  hydrogen,  and  the  latter 
refers  to  high  oxygen  concentration  regions. 

In  the  light  ion,  fast  flow  regime,  a mathematical  formulation  has  been 

worked  out  for  the  determination  of  the  pitch  Pc,  yaw  Yr  and  maqnitude 

r (4)  r 

f of  plasma  flow.  Details  are  reported  in  reference'  . The  results  are 
as  follows: 

tan  (Et  * PF)  - V cos<yt  s™£t  ~ 11  cosY0  sin  (Et  * V 

V COS  £Tt  COS  £t  - q COS  Yq  COS  (Et  + Pq) 

V sin  <5*  - q sin  Yn 
tan  Yc  1 Q 


fv2  cos2  <rt  + q2  cos2  Yq  - 2Vq  cos  <rt  cos  PQ  cos  Yr 


f2  = q2  + v2  - 2q  v (cos  cos  cos  p + sin  ^ sin  Y ) 


where 


Pn  = tan'1  )R; 
Q IT 


1*2.4  1-X  1 

|tan«c  1+XR,  . 1 
r»  L - <-»4- 


-£(t) 


Y = tan'1!  f*1  »3  cos  ^ sine"!  cos  (*+Pn+e(t))i 

tQ  jj^  , sin  £ - tan  2^  cost  ^ ) 

Rnj  ■ (Ii  - ‘j>/  (Ii  * V 

and^u  = (Rj  3 + tan  2«<tan€l  J 3 ■ tan  2#*tan€  ). 

In  these  results,  the  effects  of  spin  tilting  and  ion  sensor  skew 
orientations  have  been  taken  into  account. 

The  experimental  data  show  sharp  drop  off  behavior  in  the  ion  current 
time  variations  for  ionospheric  regions  of  heavy  ion  mass  low  speed 
plasma.  To  determine  possible  presence  of  a dominant  thermal  ion 
effect,  extra  terms  due  to  thermal  ions  were  added  to  the  denominator  of  R.  . 
while  this  effect  is  cancelled  out  in  the  numerator.  Three  simultaneous  trig- 
onometric equations  resulted  and  had  to  be  solved  simultaneously  at  every 
second.  They  are  of  the  form:  - 

2 2 

cos oi  (R2  j cos  t - sin  X)  - sin©*  sinTfcosTf (R2  j + D 

+ R-  , C = 0 

1 2 p 

cos©*,  (R2  3 cos  X - sin  X)  + sin©*  sinycosy(R2  3 + 1) 

*R2.3  C -° 

cosoi  (2R3  i cos2X  ” R3  1 ) ” 2 sin o<  sinXcosy 

+ r3>1  C = 0 

It  was  found  that  thermal  ion  effect  is  not  predominant  . 

To  investigate  the  density  and  temperature  of  the  ion  plasma  in  the 

low  speed  regime,  it  is  necessary  to  perform  nonlinear  fitting  to 

(51 

the  ion  current  equation' 


*)3/2  ( { { 

Y w a ' L 


I = Ne  A £)3/2  jf  jf  l V2exP{-a^2  + 

Vy  + + q2  - 2q  (Vg  cos  ff  + V sinO)  [ c 

dvy  dvz  1 ' 


Integrating,  we  have 


I » % Ne  A q cos  0-  £ 1 + erf  (x)  + exp  (-x^/xlT^j 
x ■ q a **  cos  0 * Uq/cif15)  cos  O’ 


where 

and  erf  (x)  Is  the  error  function  defined  as 

.x 

/o  \ / 

erf(x) 


fe)/. 


exp(-t^)  dt 


Here,  q Is  the  plasma  flow  velocity  which.  In  low  altitude  and  equatorial 
regions.  Is  mainly  due  to  the  relative  motion  of  the  satellite  In  the 
plasma  environment; cos  0 Is  the  angle  between  flow  direction  and  sensor 
normal . 


The  parameters  to  be  fitted  are  ot(*Ne  A)  and  p (=  a1*).  The  former 
(•<)  measures  the  Ion  density  and  effective  sensor  aperture  area  and 
Is  responsible  for  the  height  of  the  Ion  current  variations.  The 
latter(B)  Is  also  / m ~ which  measures  the  plasma  velocity,  or  Ion 

mass  to  Ion  temperature  ratio,  p influences  the  profile  of  the  Ion 
current  curves.  A typical  nonlinear  fitting  of  Ion  current  for  a 
satellite  spin  is  shown  In  Figure  4. 
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13.0  DMSP  Topside  Plasma  Monitor  (Phase  I) 


Initiator:  R.  Sagalyn 

Project  No.:  8617  Problem  No.:  4836 

The  Topside  Plasma  Monitor  to  be  flown  on  Defense  Meterological  Systems 
Project  (DMSP)  satellites  consists  of  two  sensing  systems;  one,  the  electron 
sensor,  measures  ambient  electron  density,  electron  temperature  and  vehicle 
potential  and  the  other  system,  the  ion  sensor,  measures  the  ambient  density 
of  each  ion  specie,  the  average  ion  temperature  and  the  vehicle  potential. 

The  observed  parameters,  measured  as  a function  of  time  along  the  vehicle 
trajectory,  are  to  be  combined  with  ephemeris  data  to  give  plasma  scale 
height  and  faF£  as  a function  of  position  and  density  irregularity  frequency 
spectra.  This  Supplementary  Sensor  Ion-Electron  (SSIE)  system  is  to  be 
used  in  cooperation  with  Global  Weather  Center,  Offutt  AFB. 

This  report  outlines  procedures  planned  as  a result  of  analyses  to  date, 
both  for  the  data  processing  system  and  for  determination  of  electron  and 
ion  parameters.  Preliminary  studies 

with  sample  data  are  anticipated  to  evaluate  the  one  and  two  specie  analy- 
sis procedures  for  H+,  He+  and  0+.  Stability  of  the  solution  technique 
using  a restricted  number  of  sweep  data  points  is  essential  to  providing  a 
reliable  on-line  analysis  package  for  determining  specie  densities, 
temperatures,  vehicle  potential,  and  eventually  the  plasma  scale  height 

and  f_F^. 

o c 

13.1  Data  Processing  System 

The  data  processing  system,  electron  and  ion  mass  and  temperature  analyses, 
and  the  file  structures  are  defined  for  operation  in  a real  time  environment 
on  a UNIVAC  1110  coi.1; at  Offutt  AFB.  Initial  development  and  follow- 
on  global  mapping  work  is  to  be  conducted  on  the  CDC  6600  system  at  AFGL. 
Compatibility  of  the  program  between  the  two  systems  will  be  maintained 
by  processing  36  bits  of  data  from  the  UNIVAC  1110  in  one  60  bit  word  on 
the  6600.  Requisite  conversion  of  floating  point  numbers  will  account 
for  the  shorter  characteristic  and  mantissa,  and  the  sign  plus  magnitude 
arithmetic  of  the  1110.  ^ ^Other  special  conversion  considerations  include 
the  use  of  IOWAIT  and  FLD  instructions  with  GWCI0,  instead  of  WRITMS,  READMS 
and  the  character  and  bit-processing  instructions  (MXGETX,  LBYTX,  etc.) 
with  the  6600  SCOPE-Fortran  System. 


Data  processing  will  consist  of  read  out  of  multiples  of  28-word  sectors 
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Univac  1110 


CDC  6600 


Store  backward 
in  time.  Jumps 
or  Overlap  to 
be  Edited. 


Contents 

1)  Directory 

2)  Read-outs  of  orbits 

Orbit  Read-out  contains 
(multiples  of  28-word  sectors) 

a)  Info  Block  - Satellite  & time 

b)  Minute-by-minute  blocks 

i)  Ephemeris  group  - 32  words 
ii)  Data  - 6 words/second  (packed) 


Same  as  1110 
except  for 
IOWAIT,  FLO. 
Floating  Pt. 
differences. 


Analysis  Results 

Nj,  N2.  T.  T£,  N,  Mj,  ^se>  Hp,  fQF2 


Similar  to  input  file 
except  minute-blocks 

i)  Ephemeris  and  Sweep 

Analysis  group  - 60  words 
ii)  Data  as  before 


Figure  1.  Defense  Meteorological  Systems 
Project  SSIE  Sensor  Processing  System 
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from  the  input  circular  file,  editing  and  upacking  of  minute-by-minute 
blocks  of  ephemeris  and  instrumentation  data,  forward  time  handling  of  the 
backward  recorded  data,  augmentation  of  the  ephemeris  group  in  each 
minute  block  by  the  results  of  the  analysis  of  the  electron  and  ion  sweeps, ^ 
and  recreation  of  an  output  circular  file*  Observed  vehicle  potential 
and  densities  from  each  system  will  be  compared  as  a check  on  satisfactory 
instrument  performance.  Debug  and  standard  output,  archiving  capability, 
and  mapping  programs  are  required.  A chart  outline  of  files  and  functions, 
applicable  either  on  the  UNIVAC  1110  or  on  the  CDC  6600,  is  presented  in 
Figure  1. 


13.2  Determination  of  Electron  and  Ion  Parameters 


The  raw  data  obtained  from  the  electron  and  ion  sensor  in  sweep  mode 

consist  of  current  as  a function  of  applied  voltage.  From  this  data  electron 

and  ion  density  and  temperature  must  be  deduced,  as  well  as  the  vehicle 

(21 

potential  with  respect  to  the  plasma'  . Variations  in  densities  between 
sweep  modes  are  obtained  by  comparison  of  current  at  a fixed  applied 
voltage  to  that  obtained  in  the  sweep  mode  of  the  appropriate  sensor  for 
the  same  voltage. 


13.2.1  Electron  Sweep  Mode 


The  electron  sensor  current  and  slope  as  a function  of  applied  voltage  are 


e 0 (Retarding! 
X =Ae  <K  N*a  gxp.__(-xa) 

2 /tT 


e 0 < O(Accelerating) 
(l  - *2) 


Ae«N  a 

e 


t/i r 


3i 

*0i 


_ -Ae  ©< 

Ty 


r ] 


-Ae  * N*a 

Tpw 


where 


2 e 0 


0 = potential  of  sensor  with  respect  to  plasma  =0{.F+0  (volts) 

_ P 

0 SE  = vehicle  potential  with  respect  to  plasma  (volts) 

= sensor  potential  with  respect  to  vehicle  (volts) 

P 2 2 

A = sensor  surface  area  = 47Tr  (m  ) where  r is  sensor  radius  (m) 
e = electronic  charge  = -1.6021  x 10‘19  (M.K.S  u's) 

oc  = sensor  transparency  =0.8 
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Ng  = ambient  electron  density  (m  ) 
mg  = electron  mass  (Kg) 

a = fekTe  = most  probable  electron  speed  (m/s) 

/ "e 

Tg  = electron  temperature  ( K) 

k = Boltzmann  constant  = 1.38054  x 10~2^  (M.K.S  u's) 

This  slope  is  expected  to  be  constant  for  accelerating  voltages.  In 
actual  practice  the  slope  will  increase  with  increasing  0 due  to 
Debye  Shielding.  Thus  it  should  reach  a minimum  at  £=  0 or  0p  = - 0^. 
The  vehicle  potential  0$E  is  therefore  obtained  by  locating  this  minimum. 
The  electron  temperature  and  density  are  obtained  from  the  logarithmic 
slope  of  the  negative  of  the  current  at  this  point  and  the  current 
respectively. 


13.2.2  Ion  Sweep  Mode 

The  ion  sensor  current  as  a function  of  applied  voltage  is 


For  e 0£  0 (Retarding) 
N, 


I = AeJk_gc  £ 
2 


i=r 


1 + erf  (x^ ) +d4 


exp 


where 


x 

2 

exp  (-z  ) d z - definition 

0 = potential  of  sensor  with  respect  to  plasma  = 0g  + 0p(  volts) 

0S  = potential  of  vehicle  with  respect  to  plasma  (volts) 

0 = sensor  potential  with  respect  to  vehicle  (volts) 

P 2 2 

A = sensor  aperture  area  = TTr  (m  ) where  r is  aperture  radius  (m) 

e = electronic  charge  = 1.60210  x 10"19  (M.K.S.  u's) 

<x.  = sensor  transparency  = 0.59049 

i.L  O 

N,  = ambient  ion  density  of  i constituent  (m”J) 

• XL. 

iHj  = mass  of  iLn  constituent 


erf  (x) 


=AJo 
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T 

k 

j 

V. 


= most  probable  speed  of  i^*1  constituent  (m/s) 


average  ion  temperature  of  all  constituents  ( K) 
Boltzmann  constant  = 1.38054  x 10’23  (M.K.S.  u's) 
number  of  ion  constituents 
vehicle  speed  (m/s)  = 7.437172213  x 103  m/s 


For  e0<O  (accelerating) 

Put  0 = 0 in  retarding  equation  whence  x^  = V$/ai 

and  I is  constant  (not  a function  of  0 ) so 

P 


Only  two  species  H+  and  0+  are  considered  since  they  are  expected  to  dominate 
at  the  altitude  of  interest  (835  km).  From  the  above  we  expect  to  find 
a minimum  in  the  slope  of  I vs  0 for  each  of  these  species  present.  This 

r 

corresponds  to  the  voltage  just  sufficient  to  repel  a stationary  ion 
(relative  energy  H m^V$  ) from  the  detector.  When  both  ions  are  present 
the  two  minima  will  be  widely  separated  and  the  slope  at  each  can  be 

assumed  to  be  due  only  to  the  one  ion.  The  procedure  to  be  followed 

depends  on  whether  one  or  two  minima  in  the  slope  are  detected. 

13.2.2.1  Two  Ion  Case 


In  this  case,  indicated  by  the  presence  of  two  minima  in  the  slope  of 
current  applied  voltage,  the  temperature  and  0+  density  are  obtained 
from  solution  of  simultaneous  equations  for  the  current  and  slope  at 
the  0+  minimum,  since  H+  contributes  to  neither  at  this  point.  This 
yields 


n2 

= 16m 3 T*  r.'< 

. rvs2  . 2 

I2e 

Ae2o<  1.2 

L 4 16mpS2  J 

T 

= iA  e^n/  e 

2 

2TTk  16.mpS2z 

is 

the  current, 

the  negative 

slope 
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The  vehicle  potential  is  obtained  as 


the  proton  mass. 


#s 


1 

? 


F 


mV 

— P—L 


2e 
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where  and  0^  are  the  locations  of  the  two  minima.  The  density  of 

H+  is  then  directly  obtainable  from  the  current  at  0 = - 0 . The 

p s 

average  ion  mass  M may  also  be  obtained. 

13.2.2.1  One  Ion  Case 

In  this  case  a rigorous  algebraic  solution  of  the  equations  is  impossible. 
Two  approximate  methods  are  used  to  estimate  the  ion  density 

1)  Nj  = electron  density 

2)  N,  obtained  from  current  I at  0 = 0 which  is  in  the 

1 ° p d 

accelerating  region  (I  = constant,  7f\.liP*  °),by  setting 
the  expression  in  the  brackets  = 2. 

( Xx>>  Q,  erf  (X1)-  1,  exp  (-X^)-  0) 


The  equations  for  current  and  slope  at  the  slope  minimum  (-S-|)  are  solved 
for  mass  and  temperature  yielding 


■j  = (fl  e 


T = 


(Ae  <*  N-,)2  e2 
2TTk  m]  S]2 


The  vehicle  potential  is 


13.2.2.1.1  Second  Method  of  Estimating  Mass 

Small  errors  in  the  current  I,  or  density  N,  can  lead  to  large  errors  in  the 
mass  as  obtained  above  because  of  the  subtraction  which  occurs  in  the 
denominator.  This  is  particularly  dangerous  if  the  ion  happens  to  be  H+ 
at  low  temperatures,  for  which  the  slope  minimum  will  be  steep.  Hence 
a second  method  is  proposed  based  on  the  fact  that 

^•<T  <Te 

where  Te  is  the  electron  temperature.  The  ion  present  (H+or  0+)  can  then 
be  assumed  to  be  the  one  which  yields  a temperature  closest  to  Te. 

13.2.3  Plasma  Scale  Height  and  fQF^ 

We  have  sufficient  information  from  the  M2  EL  and  M2  ION  modes  to  compute 

the  plasma  scale  height  H at  a time  centered  on  the  two  modes  which  occur 

P 


17  secs  apart,  that  Is  H„  can  be  deemed  to  be  determined  at  (t  +8.5) 

p o 

where  tQ  Is  center  of  electron  sweep,  and  re-evaluated  at  an  Interval 
of  128  secs. 

Hp  + 7*5)  s k (Te  + T)/M  mp  g (m) 

Where  g ■ acceleration  due  to  gravity  at  835  Km  altitude  (m/sec2) 

g - 7.675373268  m/sec2 

M ■ average  Ion  mass  In  AMU's 

Substitute  In  constants 

(t_  + 8.5)  = 1.083302743  (Te  + T)  Km 
p o R 

fQF2  will  be  determined  assuming  a diffusive  equilibrium  model  with  a 
specified  peak  height  of  the  F-reglon  as  a function  of  local  time, 
and  Is  given  by  the  formula 

nfi  ■ 1.24  x 104f2. 
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14.0  High  Latitude  Scintillations  Studs 


Initiator:  J.  Aarons 


Project  No.  4643 

14.0  High  Latitude  Scintillations  Studs 


Problem  No.:  4785,  4844 


In  order  to  provide  a realistic  statistical  and  dynamic  model  of  scintillation 
occurrence  at  sub-auroral  and  auroral  latitudes  a model  development  was  under- 
taken. The  primary  data  consisted  of  several  years'  measurements  at  three 
stations  of  15-minute  samples  of  scintillation  excursions  in  dB  of  the  137 
MHz  beacon  from  the  ATS-3  synchronous  satellite. 

The  initial  model  was  to  provide  mean  scintillation  excursions  as  a function 
of  month,  time  of  day,  magnetic  index,  and  solar  flux.  In  addition,  the  data 
base  would  be  used  for  statistical  distribution,  regression,  and  short  term 
forecasting  studies. 

14.1  Data  Base  Creation 

Measurement  data  was  provided  on  punched  cards  with  16  15-minute  samples  or  4 
hours  per  card.  The  scintillation  data  card  format  is  given  below: 


Column  # 1- t 

7 

10  12-59 

68-69 

70-72 

73-74  79-80 

312 

11  2X, 

IX,  1613  9X 

, A2 

13 

12  4X  , A2 

YR.MO.DA 

Period 

1 4 15min 

Sat 

Freq. 

Stat.  * DB ' 

of  day 

9 readings 

ID 

ID 

1-6 

or  0 (x  4 hrs. 

) e.g. 

A3  e.g. 

137 

1=001 

low  (whole  period) 

9=-99 

no  data  ( " " ) 

0=normal  (read  each  column) 

me  se 
CDC  I 


The  Kp  values  were  provided  by  Dr.  P.  Fougere  of  the  Space  Physics 
Laboratory  and  the  solar  flux  values  were  provided  by  Mr.  W.  Barron  (LI).  The 
ATS-3  satellite  ephemeris  was  generated  using  Spadats  element  sets  in  the 
LOKANGL  program  provided  by  the  Analysis  and  Simulation  Branch.  Sub-ionospheric 
latitude  and  longitude  at  350  Km  were  obtained  using  routine  SILL2  provided  by 
the  Ionospheric  Physics  Laboratory.  The  corresponding  geomagnetic  coordinates 
and  geomagnetic  time  were  then  obtained  based  on  G.  Gustafsson's  1970  Revised 
Corrected  Geomagnetic  Coordinate  System,  using  routines  C0RRGM2  and  MAGTIM 
developed  previously  for  the  Ionospheric  Physics  Laboratory.^ 

The  average  geographic,  geomagentic  and  geometric  parameters  for  the  sub-iono- 
spheric points  for  the  three  stations  are: 


Geographic 

Inv. 

Lat.  Long. 

Lat. 

El. 

A z. 

Narssarssuaq 

5.42  51.0 

63. 2* 

18.0° 

208* 

Goose  Bay 

48.3  61.7 

60.3° 

28.8s' 

191* 

Sagamore  Hill 

39.3  70.6 

53.5° 

40.9° 

178° 

A full  data  base  tape  comprising  each  15-minute  sample  and  environmental  parameters 
was  generated  for  each  station;  the  record  format  for  each  sample  is  described  in 
Table  1. 
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VARIABLE 


10  (station-data  source) 

SAT  ident. 

FREQ 

IY  Year 
IM  Month 
I DAY  Day 
UT  Hr. 

SI  Scintillation  Index 
"OB" 

Kp  Geomagnetic  Index  (3-hourly) 

SF27  2.7  GHz  solar  flux 

SF50  5.0  GHz  solar  flux 

SLAT  sub-ionos.  latitude 

ELONG  sub-ionos.  long.  (E  positive) 

CGLAT  sub-ionos.  mag.  lat. 

CGMT  mag.  time  (hr.) 


FORMAT 


TABLE  1.  Scintillations  Data  Base  Tape  Format 


T 


The  time  span  and  size  of  the  data  bases  are  as  follows: 


Narssarssuaq 
Goose  Bay 
Sagamore  Hill 


9/17/68  - 9/1/74 
1/1/72  - 12/31/74 
12/1/69  - 11/30/74 


146,700  + samples 
71,000  + samples 
148,000  + samples 


Goose  Bay  data  for  1974  cover  mainly  November  and  December. 


14.2  Modeling 

Analyses  were  conducted  separately  for  each  station.  In  order  to  uniformize  the 
data  for  modeling  studies  the  data  were  partitioned  into: 

12  months,  7 Kp,  3 Solar  Flux  (2.7  GHz),  and  24  UT  ranges.  The  date, 

Kp,  Sf  and  SI  readings  were  averaged  in  each  block.  A compact  file 

was  thus  made  available  for  high  speed  iterative  modeling  studies.  The 
seven  Kp  ranges  are  0-1,  1+  to  2,  2+  to  3,  3+  to  4,  4+  to  5,  5+  to  6, 

6+  and  up.  The  three  Sf  ranges  are  0 to  95,  96  to  120,  121  and  up. 

Tables  of  the  averaged  SI  were  provided  for  each  of  the  stations. 


Out  of  a maximum  possible  6048  blocks  (12  x 7 x 3 x 24),  the  averaged  files 
comprised  the  following: 

Narssarssuaq  4985  blocks 

Goose  Bay  4217  blocks 

Sagamore  Hill  5065  blocks 


The  empty  blocks  generally  correspond  to  the  highest  two  Kp  ranges,  i.e.,  5+  and 
up,  and  occasionally  the  highest  Sf  range.  Table  2 shows  the  averaged  SI  by 
hour,  solar  flux,  and  Kp  range  that  were  obtained  for  January  and  February  for 
Goose  Bay.  ’ R'  implies  absence  of  data. 
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An  extensive  search  was  conducted  to  derive  empirical  models  of  SI  for  each  of 
the  three  stations.  Analytical  forms  of  the  model  were  preferred  to  ensure 
smooth  transitions  as  a function  of  the  known  driving  factors,  viz.  day  of  year, 

Kp,  5^  and  universal  time.  These  forms  also  permitted  use  of  regression  techniques 

for  least  squares  fitting  to  the  averaged  data  file.  Program  MODEL  was  developed 

(2) 

by  adapting  the  IBM  multiple  linear  regression  program  REGRE'  . The  square  of 
the  residuals  of  all  available  averaged  data  blocks  is  minimized  when  the  multi- 
ple correlation  of  the  partial s of  all  the  coefficients  that  are  to  be  determined 
vanishes.  Selective  "fixing"  of  coefficients  is  permitted  if  the  solution  exhibits 
instability,  and  an  automated  procedure  is  available  for  rapidly  converging  to 
the  optimum. 

In  the  course  of  the  search  for  improved  fits  special  characteristics  of  the  data 
were  noted  which  suggested  elaborations  of  the  model  form.  Examples  are  the 
delayed  peak  in  the  diurnal  SI  variation  with  higher  Kp,  the  seasonal  effect  on 

diurnal  variation  amplitude  relative  to  the  average  SI,  the  seasonal  effect  on 
influence  of  Kp  and  Sf,  and  the  need  for  higher  harmonics  to  represent  the 

diurnal  variation.  Figure  1 presents  the  optimized  model  that  was  obtained  for 
Goose  Bay.  An  identical  model  form  with  different  coefficients  was  also  obtained 
for  the  other  stations^.  The  cosine  terms  for  the  seasonal  variation  omit 
21^365  and  the  cosine  terms  for  the  diurnal  variation  omit  21^24  for  convenience. 

A$  is  0.01  times  the  2.7  GHz  solar  flux;  and  time  is  in  hours  UT. 

Full  comparative  plots  of  the  averaged  data  and  the  model  were  provided  for  each  station. 
The  model  predictions  used  the  actual  averaged  data,  Kp  and  for  each  hour,  and 

are  therefore  absent  when  data  are  absent.  These  best  fit  models  may  occasionally 
predict  small  negative  Si's;  these  should  in  practice  be  set  to  a selected  minimum. 

Figures  2-5  illustrate  the  behavior  of  the  averaged  data  and  of  the  best  fitting 
model  equations.  The  features  of  increasing  scintillations  with  increasing  Kp,  and 
generally  with  increasing  solar  flux  are  evident.  The  diurnal  and  seasonal 
variations  are  also  modeled.  Standard  deviations  of  the  order  of  2 dB  are 
obtained  with  these  models. 

Subsequently,  effort  was  directed  to  bridging  the  models  for  the  three  stations 
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Figure  7.  Model  of  Average  SI  for  Goose  Bay 
Optimization  is  according  to  Table 


GOOSE  BAY  MONTH  3 

SOLAR  FLUX  0-95  96-120  121- UP 


Figure  2.  Average  SI  Data  and  Model  Plots  for  Goose  Bay  for  March 


by  going  to  a local  time  base  and  by  Including  the  geomagnetic  latitudinal  de- 
pendence. The  three  separate  models  provided  the  "data"  and  the  new  comprehensive 
model  was  developed  using  program  MODEL.  The  global  model  equation  for  the 
three  stations  combined  is  presented  in  Figure  6.  Figure  7 shows  the  plot  for 
March,  in  local  time,  for  the  cases  of  low  and  high  solar  flux  (80  and  140)  under 
four  selected  magnetic  canditions  (Kp  - .5,  2.5,  4.5  and  6.5).  Standard  deviation 
of  the  residuals  between  the  global  and  separate  models  was  less  than  0.8  dB, 
indicating  only  a moderate  further  deviation  from  the  original  sample 
data.  Greater  detail  on  various  aspects  of  this  study  is  presented  in  Reference  3. 


14.3  Distribution  of  Scintillation  Index  Readings 

Distribution  analyses  were  conducted  for  each  of  the  three  stations.  The  effect 
of  the  solar  flux  was  not  evaluated,  but  partitioning  by  month,  Kp  and  hour  was 
retained.  Within  each  partition  the  average  SI  was  calculated,  and  the  percentage 
of  readings  falling  in  consecutive  ranges  was  determined.  Significant  features  of 
the  distributions  have  been  obtained  by  examination  of  the  tabulations,  and  are 
discussed  below.  For  selected  average  SI,  the  percentage  SI  readings  found  in 
the  ranges  0-1,  1-3,  3-6,  6-9,  9-12  and  12  and  up  are  presented  in  Table  3. 

For  each  station,  no  marked  variations  in  the  distributions  are  evident  for 
different  months,  Kp  ranges  or  time  of  day.  Thus,  given  an  average  SI 

value  for  any  time-Kp  partition  for  a specific  station,  the  expectation  of  any 

range  of  readings  may  be  predicted.  Differences  in  distribution  patterns  among 
the  three  stations  is  primarily  due  to  the  consistently  higher  SI  readings  at 
Narssarssuaq  on  the  one  hand,  and  the  extremely  low  activity  at  Sagamore  Hill  on 
the  other. 


Averaged  vs.  Median  Data 

The  actual  data  plotted  in  Figures  2 through  5 upon  which  the  model  is  based,  is 
the  averaged  scintillation  data  in  dB.  It  was  desired  to  compare  this  mean  with 
the  calculated  median  for  the  same  time,  Kp,  and  solar  flux  blocks.  For  this 

purpose  the  full  data  base  was  sorted  by  scintillation  index  within  each  block, 
and  a program  SIMED  obtained  the  medians  for  creation  of  a packed  file  similar 
to  the  averaged  file.  Figure  8 illustrates  both  the  mean  and  median  dB  values 
for  Goose  Bay.  This  data  Is  for  the  same  time  period  (March)  as  described  in  the 
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Figure  6.  Global  High  Latitude  Model  of  Average  SI  Based  on  ATS-3  Data  From 
3 Stations 


„ . TIME-UT 

e 8.  Average  and  Median  SI  for  Goose  Bay  for  March 


earlier  figures.  It  was  seen  that  no  substantial  differences  occur  between  the 
averaged  and  median  dB  values  at  any  of  the  stations,  except  that  for  low 
scintillation  activity  median  values  tend  to  be  lower  than  the  averaged  values. 

i 

14.4  Other  Analyses 

Some  intermediate  analyses  were  conducted  during  the  course  of  this  study.  For 
each  station,  month,  and  hour  the  averaged  file  was  used  to  run  a multiple 
linear  regression  for  scintillation  index  as  a function  of  Kp  and  2.7  GHz  solar 

flux  (FL2).  The  coefficients,  correlations,  and  the  multiple  correlation  was 
output  for  each  Interval.  Table  4 is  a sample  of  the  output  by  hour  for  the  first 
six  months  for  Goose  Bay.  Plots  of  the  correlation  coefficients  were  also  gener- 
ated. In  general  the  correlation  with  Kp  is  much  greater  than  with  solar  flux. 
Multiple  correlations  are  consistently  high  (about  0.8)  for  Narssarssuaq,  around 
0.65  for  Goose  Bay,  and  about  0.5  for  Sagamore  Hill. 

Separate  regression  analyses  also  showed  greater  correlation  with  Kp  during 
periods  of  higher  magnetic  activity. 

A forecasting  analysis  attempted  to  predict  early  1975  Goose  Bay  scintillations 
that  were  not  included  in  the  data  base.  The  results  showed  unexpectedly  high 
SI  observations,  possibly  because  of  new  instrumentation  sensitive  to  high  SI 
excursions. 
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15.0  Real  Time  Signal  Processing  With  CSP-30 


Initiator:  J.  Schindler 

Project  No.:  5635  Problem  No.:  4794 

15.1  Introduction 

The  purpose  of  this  project  is  to  develop  assembly  language  real-time  software  for 
the  CSP-30  digital  signal  processor.  The  software  will  be  used  in  a simulation 
experiment  designed  to  test  recently-developed  filtering  techniques^ Vor  maximi- 
zing si gnal-to-cl utter  ratio  when  detecting  a moving  target  with  moving  radar. 

The  processor  will  receive  16  bit  complex  signals  from  each  of  several  antenna 
elements  (eventually  8)  at  a pulse  repetition  rate  of  1.8  KHz.  For  each  antenna 
element  and  pulse  a signal  is  received  from  each  of  up  to  16  contiguous  range 
bins.  For  each  range  bln  and  pulse  a linear  combination  of  the  signals  for 
the  different  antenna  elements  are  formed  which  define  the  radar  patterns. 

These  are  then  filtered  and  the  outputs  displayed  on  a storage  display  scope 
along  with  other  pertinent  information. 

15.2  Maximum  Signal  to  Clutter  Processing 

The  Doppler  shift  of  a signal  received  by  a moving  radar  (velocity  vQ)  from  a 

stationary  clutter  element  at  angle  e with  respect  to  the  velocity  vector  is(^ 

VS**  sin  0 (1) 

Where  X Is  the  wavelength  of  the  signal.  Thus  for  equal  transmitted  and  received 
antenna  patterns  P (e)  the  doppler  clutter  spectrum  is 
D Ud)  - P2  [e(wd)j  * P2  ^sin'1  ^dj 

The  presence  of  a moving  target  within  a given  Interval  of  the  doppler  spectrum 
may  be  detected  by  constructing  a received  antenna  pattern  with  a null  over  the 
appropriate  angular  range  to  exclude  ground  clutter  contributing  to  this 
interval  and  filter  the  resulting  signal  to  pass  frequencies  within  this  Interval. 
Reference  1 describes  a method  of  selecting  the  coefficients  defining  either 
the  antenna  pattern  or  filter  to  maximize  the  signal  to  clutter  ratio  by  deriving 
an  expression  for  It  which  is  quadratic  In  the  set  of  coefficients  (antenna  or 
filter)  which  are  to  be  varied.  The  results  are  as  follows: 
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(1)  Radar  coefficients 
' , \ 


X = 


’1 


X J 
N 


where  is  the  complex  coefficient  of  the  n^h  element  in  the  radar  pattern 


X = B"1  a 


where 


N N 

L I 


Bnp  n-t.  mn 


F[  Pt  (sin'y  ) f(sin-^)  (1-jft* 


J dy  exp  £-j  (m-q)  dkyj 


sqp  / dy  exp 


1 = $+i* 


wi  th , 


S = S + S ' 
mn  mn  mn 


Smn  = field  scattering  coefficient  from  mth  radar  element  to  nth  radar 
element 
k = 2 IT/  a 

d=  spacing  between  antenna  elements 

f(6)  = cowon  pattern  rtf  individual  antenna  elements 

* = complex  conjugate  operation 

t = complex  conjugate  transpose  operation 


3 


202 


F(f)  = doppler  filter  frequency  response 
Pt  (0)  = transmitted  antenna  power  pattern 
= expfejxkdsin©t) 

0t  = angle  of  target 


(2)  Filter  Coefficients  <*.,• 
-1_ 

ot  =1  a, 


Where 


B„  = 
np 


j Pr (9)Pt (9)  exp  ^2jkv0(n-p)Tsine]  de 


an  * exp  (-jt^nT) 

With 

pr(9)  = received  antenna  pattern 
T = time  spacing  between  samples. 
ut  = doppler  frequency  of  target 


15.3  Hardware 

Figure  1 shows  the  probable  configuration  in  which  the  CSP-30  will  be  used  in 
this  experiment.  The  CSP-  3&»1  s a 16  bit  digital  computer  with  signed  fixed 
point  arithmetic  and  a basic  cycle  time  of  100  nanoseconds.  It  contains  4096 
words  of  high  speed  integrated  circuit  (IC)  memory  and  32767  of  low-speed 
(600  n sec)  core  memory.  In  addition  there  exists  an  accumulator  file  of  32 
registers  and  a super  accumulator  or  A register.  Background  instructions  are 
available  to  initiate  transfer  of  data  between  core  and  the  higher  16  registers, 


partially  concurrent  with  the  execution  of  subsequent  instructions. 

Programmed  one  word  data  transfers  may  be  Initiated  to  and  from  the  tapes  and 
teletype.  The  digital  radar  signals  are  to  be  input  as  blocks  through  the  core 
direct  memory  access  channel  (CDMA).  The  storage  display  scope  contains  a 
1024x1024  point  screen  on  which  display  may  be  stored  at  the  option  of  the 
program.  Hardware  facilities  include  plots  of  points,  cursors,  x and  y axes. 

Priority  interrupt  logic  and  status  words  (not  shown)  permit  overlapping  of  I/O 
and  other  operations.  CDMA  operations  may  be  overlapped  with  any  other  operation 
not  accessing  core  memory. 


' > 

A direct  communication  between  the  operator  and  computer  is  possible  through  the 
master  switch  register  (MSR),  a set  of  16  console  switches,  numbered  0 through  15. 

A special  high-speed  (MAP)  processor  to  be  provided  by  CSPI  will  facilitate 
the  real-time  processing  to  be  done. 


15.3.1  Timing  Considerations 

Table  1 shows  significant  timing  data  which  will  effect  this  experiment.  The 
complex  multiplication  operations  include  4 single  precision  multiplications, 

2 addition  and  any  memory  access  required  (other  than  core).  The  CDMA  data 
transfer  time  is  for  the  case  in  which  the  CDMA  has  exclusive  access  to  core 
memory  (an  option  available  to  the  program).  Timing  information  for  data  trans- 
fers to  and  from  the  MAP  processor  is  not  available  as  of  this  writing. 

It  is  expected  that  several  processes  might  be  taking  place  simultaneously  to 
make  optimum  use  of  the  facility,  such  as, 

A.  Process  the  data  (form  radar  patterns  and  filter)  for  one  pulse 

B.  Input  data  for  next  pulse 

C.  Output  results  from  previous  pulse 

Based  on  Table  1 the  following  time  estimates  are  obtained: 

A.  Form  p antenna  patterns  for  m range  bins  n antenna  elements  and  filter 
with  f point  filter 

tA  - mp(n+f)tM 

B.  Input  next  pulse  to  IC  memory 

tj  ■ 4.4mn/«  sec 

C.  Display  results  of  previous  pulse  (one  point  for  each  antenna  pattern  of 
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each  range  bin) 

tQ  * 30  mp  y.  sec 

The  advantage  of  using  the  MAP  processor  for  (A)  is  evident.  Even  if  this  is 
the  case,  however,  it  is  apparent  that  at  most  3 radar  patterns  can  be  formed 
for  the  planned  8 antenna  elements  and  16  range  bins,  in  the  500  fl  sec  pulse 
interval.  Thus  the  number  of  range  bins  may  have  to  be  reduced  to  obtain 
better  Doppler  resolution.  The  timing  requirement  tp  could  be  reduced  by  dis 

playing  only  the  largest  outputs. 


TABLE  1.  Approximate  CSP-30  Timing  Data 


complex  Multiplication,  tM 
(CSP-30) 


7 p sec 


complex  Multiplication,  tM  1 ji  sec 

(MAP) 


CDMA  Data  Transfer 


1 y-  sec/word 


Core  to  accumulators le  .6  >»  sec/word 

Data  Transfer 


Accumulator  File  to  .6  yu  sec/word 

IC  Data  Transfer 


Display 


30  }i  sec/point 
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15.4  Fast  Fourier  Transform  Display 


To  gain  familiarity  with  the  computer  a program,  described  in  Figures  2-5,  was 
written  to  compute  and  display  the  discrete  Fourier  transform  of  some  simple 
functions.  The  discrete  Fourier  transform  of  an  N point  time  sequence  f.  is 


k * 


-2TTj(i-l  ) 


To  compute  F^  the  program  makes  use  of  the  Fast  Fourier  Transform  (FFT)  algorithms 
Incorporated  Into  a signal  processing  software  supplied  with  the  computer.  The 
"VIDEO"  function  referred  to  in  Figures  2 and  4 is  described  in  Reference  4. 
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Figure  2.  Fast  Fourier  Transform  Display: 

Description  of  Options.  Here  j=  V-l. 

I is  point  number  (from  1 to  1024), 

N is  in  units  of  PRF/1024,  and  the  pulse  repetition  (sampling) 
frequency  PRF  is  arbitrary. 
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Figure  3.  Console  Master  Switch  Register  (MSR)  Settings  to  Implement  Various  Options 
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Figure  4. 
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16  Point  Complex  Simulation  Signal  Generated  by  Subroutine 
VIDEO  (Reference  4) 
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